NUMERICAL STUDY OF THERMOVISCOUS EFFECTS IN . . . PHYSICAL REVIEW E 90, 043016 (2014)
FIG. 1. (Color online) (a) Sketch of a physical acoustophoresis
setup with a silicon microfluidic chip on top of a piezoactuator such
as in Ref. 30. (b) Sketch of the model system for the numerical
scheme with the viscous fluid domain surrounded by hard walls. The
thick arrows indicate in-phase oscillating displacement of the left and
right walls. By default, we set the width w = 380 μm and the height
h = 160 μm as in Ref. 30.
calculated as
ρ1(r,t ) v1(r,t ) = 12
Reρ1(r,0) v
∗
1 (r,0), (17)
where the asterisk denote complex conjugation.
III. NUMERICAL MODEL
The numerical scheme solves the governing equations for
the acoustic field inside the two-dimensional water domain
of a rectangular microchannel cross section, whereas the
vibrations in the surrounding chip material and piezotransducer
are not modeled. The water domain is surrounded
by immovable hard walls, and the acoustic field is excited
by oscillating velocity boundary conditions, representing an
oscillating nm-sized displacement of the walls. A sketch of
the physical system and the numerical model is shown in
Fig. 1.
A. Governing equations
The governing equations are solved using the commercial
software COMSOL MULTIPHYSICS 31. To achieve
greater flexibility, the equations are implemented through
mathematics-weak-form-PDE modules and not through the
built-in modules for acoustics and fluid mechanics. In contrast
to our previous work 20, the second-order Eqs. (13), (14a),
and (15c) are implemented in the flux density formulation,
which by partial integration avoids the less accurate secondorder
derivatives appearing in the body-force formulation.
To fix the numerical solution for the second-order mass and
momentum conservation equations, the spatial average of the
second-order pressure is forced to be zero by a Lagrange
multiplier.
B. Boundary conditions
The first-order acoustic fields are solved in the frequency
domain for a driven system, in which energy is added to
the system by an oscillating velocity boundary condition and
lost by thermal conduction through the walls. The walls are
modeled as hard thermal conductors with infinite acoustic
impedance and infinite thermal diffusivity. This approximation
is reasonable given the parameter values listed in Table II. In
the numerical model, this is implemented by zero velocity and
TABLE II. Acoustic impedance and thermal diffusivity for water,
silicon, and Pyrex glass at room-temperature values from Ref. 29.
Acoustic impedance Thermal diffusivity
Material (106 kg m−2s−1) (10−7 m2s−1)
Water 1.5 1.4
Silicon 20 920
Pyrex glass 17 6.3
ambient temperature at the walls,
T = T0 on all walls, (18a)
v = 0 on all walls, (18b)
n · v1
= vbc(y,z)e
−iωt added to actuated walls, (18c)
n · v2
ρ1
ρ0
= −
added to actuated walls. (18d)
(n · v1)
Here, n is the outward pointing surface normal, and
Eq. (18d) ensures zero mass flux across the boundary. With
the fixed temperature at the channel boundaries, we focus on
the thermal properties of the acoustic effects, and we do not
model any effects of heat generated by the piezotransducer.
In the related experiments, the chip and the piezotransducer
are thermally stabilized by a Peltier cooling element beneath
the piezotransducer, which keeps the chip at a constant
temperature using a temperature sensor and a PID controller;
see Ref. 30 for a detailed explanation.
It is not trivial how to apply the oscillating velocity
boundary condition. In our model, we wish to excite the
horizontal half-wavelength resonance, which at the top and
bottom walls leads to viscous boundary layers and the
generation of streaming flow. To avoid direct influence on this
flow from the actuation, we therefore choose to actuate only
the left and right walls at y = ±w/2. Moreover, an optimal
coupling to the half-wavelength resonance is obtained by
choosing the proper symmetry of the actuation, and therefore
in terms of the components vy1 and vz1, the boundary condition
on v1 becomes
vy1
±w
2
,z
= vbce
−iωt, vz1
±w
2
,z
= 0, (19)
where vbc
= ωd is the amplitude of the actuation in terms of
the displacement d, with d = 0.1 nm in all simulations. This
velocity boundary condition is well defined and yields results
consistent with experiments 14.
C. Convergence analysis
The weak form equations along with the boundary conditions
are solved on a two-dimensional triangular mesh using
the finite-element method; see Fig. 2. The resolution of the
physical field is determined by the spatial resolution of the
mesh and the polynomial order of the basis functions used to
represent the field in each node in the mesh. To test the validity
of the numerical model, we first check that the numerical
solution has converged, i.e., ensuring that further refining of
the mesh does not change the solution significantly.
043016-5