Simulating acoustically-actuated flows in complex microchannels using the volume penalization technique

TL;DR

This work develops a perturbation-based volume-penalized solver for acoustofluidic flows in geometrically complex microchannels. By decoupling the problem into a time-harmonic first-order system and a time-averaged second-order system, and enforcing boundary conditions via a zero structure velocity and a Stokes-drift–driven second-order boundary, the method solves the first-order equations with a direct Helmholtz solver and the second-order equations with a projection-based preconditioner within an FGMRES loop. Contour-integral techniques compute the acoustic radiation force, and extensive 2D test cases demonstrate excellent agreement with body-fitted grid results, establishing suitable penalty factors and smeared-interface widths for accurate representation of boundary conditions and gradients. The approach offers scalable performance on large grids and complex geometries, providing a path toward efficient, immersed-boundary acoustofluidic simulations and enabling exploration of streaming phenomena without re-meshing. The work also outlines limitations and future directions toward 3D implementations and moving boundaries, highlighting ongoing challenges in scalable solvers for coupled Helmholtz equations.

Abstract

We present a volume penalization technique for simulating acoustically-actuated flows in geometrically complex microchannels. Using a perturbation approach, the nonlinear response of an acoustically-actuated compressible Newtonian fluid moving over obstacles or flowing in a geometrically complex domain is segregated into two sub-problems: a harmonic first-order problem and a time-averaged second-order problem, where the latter utilizes forcing terms and boundary conditions arising from the first-order solution. This segregation results in two distinct volume penalized systems of equations. The no-slip boundary condition at the fluid-solid interface is enforced by prescribing a zero structure velocity for the first-order problem, while spatially varying Stokes drift -- which depends on the gradient of the first-order solution -- is prescribed as the structure velocity for the second-order problem. The harmonic first-order system is solved via MUMPS direct solver, whereas the steady state second-order system is solved iteratively using a novel projection method-based preconditioner. The preconditioned iterative solver for the second-order system is demonstrated to be highly effective and scalable with respect to increasing penalty force and grid resolution, respectively. A novel contour integration technique to evaluate the acoustic radiation force on an immersed object is also proposed. Through test cases featuring representative microfluidic geometries, we demonstrate excellent agreement between the volume penalized and body-fitted grid. We also identify suitable penalty factors and interfacial smearing widths to accurately capture the first- and second-order solutions. These results provide first-of-its-kind empirical evidence of the efficacy of the volume penalization method for simulating acoustic streaming problems that have so far been analyzed using body-fitted methods in literature.

Figures (13)

• Figure 1: Schematic of a 2D staggered Cartesian grid. The left half shows the coordinate system for the staggered grid with $\Delta x$ and $\Delta y$ being the grid spacing along the $x$ and $y$ directions, respectively. The right half illustrates a single grid cell with (first- and second-order) velocity components $u$ ($\blacksquare$) and $v$ ($\blacktriangle$) approximated at the cell faces and scalar pressure approximated at the cell center ($\bullet$).
• Figure 2: Schematic representation of contour selection for computing $\bm{\mathrm{F}}^\textrm{rad}$ on the immersed object, denoted in orange. First, an arbitrary smooth contour $S_\textrm{int}$, not necessarily aligned with the Eulerian grid, is considered. Subsequently, the integration in Eq. \ref{['eq:frad_final']} is performed by representing $S_\textrm{int}$ as a stair step contour $S_\infty$ aligned with the grid cell boundaries via a level set field $\psi$.
• Figure 3: Flowchart summarizing the solution algorithm.
• Figure 4: Schematic of a fluid-filled, rectangular microchannel comprising a rigid cylinder, denoted in orange. The channel dimensions are taken as $W = 150µm$ and $H= 40µm$. The immersed domain is circular with radius $R_0 = 10µm$, and its center is located at $(W/4, H/2)$. The device is acoustically-actuated by prescribing a rectilinear actuation at the left and right walls ($\Gamma^\textrm{w} \cup \Gamma^\textrm{e}$), while the top and bottom walls ($\Gamma^\textrm{n} \cup \Gamma^\textrm{s}$) are taken to be fixed.
• Figure 5: Convergence analysis for the acoustic radiation force $\bm{\mathrm{F}}^\textrm{rad}$ experienced by a rigid cylinder immersed within a rectangular acoustofluidic device (see the schematic in Fig. \ref{['fig:geom_circle']}) Each panel shows the variation of error in $\bm{\mathrm{F}}^\textrm{rad}$, where the error is assessed by comparing against the body-fitted grid solution. (a) plots the error with respect to grid refinement $h$, (b) plots the error with respect to penalty factor $p_\kappa$ and (c) plots the error with respect to the number of smearing cells $n_{\textrm{cells}}$. The tables containing the radiation force data for each case corresponding to Fig. \ref{['fig:convg']}(a-c) are provided in \ref{['appendix:table-for-arf-data']}.