Numerical Determination of the Secondary Acoustic Radiation Force on a Small Sphere in a Plane Standing Wave Field

Two numerical methods based on the Finite Element Method are presented for calculating the secondary acoustic radiation force between interacting spherical particles. The first model only considers the acoustic waves scattering off a single particle, while the second model includes re-scattering effects between the two interacting spheres. The 2D axisymmetric simplified model combines the Gor’kov potential approach with acoustic simulations to find the interacting forces between two small compressible spheres in an inviscid fluid. The second model is based on 3D simulations of the acoustic field and uses the tensor integral method for direct calculation of the force. The results obtained by both models are compared with analytical equations, showing good agreement between them. The 2D and 3D models take, respectively, seconds and tens of seconds to achieve a convergence error of less than 1%. In comparison with previous models, the numerical methods presented herein can be easily implemented in commercial Finite Element software packages, where surface integrals are available, making it a suitable tool for investigating interparticle forces in acoustic manipulation devices.


Primary radiation force potential
To obtain Equation (8)

Secondary radiation force potential
The derivation of this section follows the steps of Silva and Bruus.
According to Gorkov's potential theory, the acoustic radiation potential of any arbitrary field, except a plane travelling wave, can be obtained as: Moreover, in our case, the total velocity potential is the sum of the velocity potential of the external field and the rescattered field: Here the first term corresponds to the primary radiation potential, the last term is small compared to the second and third terms. Moreover, For the second term of Equation (S8), first ∇ ⋅ ( ) = cos cos ℎ + − = cos cos ℎ − (S12a) such that the second term of Equation (S8) is , | = , cos cos ℎ − (S12b)

Its real part being
Re , = − , cos cos ℎ + (S12c) As the first term of the secondary radiation potential, Equation (S5b), depends on the velocity potential of the external field, which is real, and the real parts of the scattered velocity potential are given by Equations (S11b) and (S12c) The second term of the secondary potential, Equation (S5b), can be calculated by splitting the scattered potential: Equation (S14b), requires the calculation of some auxiliary terms (the gradient of the scattered velocity potential): S15b) and consequently, the first term of Equation (S14b) after simplification and using 2 cos − sin = (S17c) Equation (S17b) can be written as and now Equations (S5b), (S13) and (S18) yield the secondary radiation potential: Note that this force is directly proportional to the monopole scattering coefficient of the scatterer particle.

Secondary radiation force in the polar direction
The derivatives of the different terms containing : and on substitution to Equation (3b) and Equation (10) As all terms contain sin or sin 2 , the above force goes to zero when = 0 However, when = /2, only terms cos , sin 2 or 1 + cos 2 disappear, leaving which is only zero at either the nodes or antinodes, where sin 2 ℎ = 0.

Mesh convergence analysis
Mesh convergence analysis was carried out using a uniform and a non-uniform meshing to assess the convergence speed of the numerical method and verify its robustness. In both cases, we use a scaling parameter mesh_size to define the minimum and maximum discretization steps. For the uniform mesh, the maximum mesh size is given as /mesh_size , while the minimum mesh size as /(2 ⋅ mesh_size). For the non-uniform mesh, the scattering and probe particles are meshed using the above minimum and maximum values, but the other domains (the fluid domain and the PML) are meshed using a coarser mesh, with minimum and maximum size of /(0.6 ⋅ mesh_size) and /(1.2 ⋅ mesh_size), respectively.
As the polystyrene particle in water has a lower contrast compared to the polystyrene particle in air, and consequently its relative secondary radiation force is lower, we chose this former case to analyze the mesh convergence. The results for various cases (node, antinode) and along different directions (radial and z) can be seen in the following Figs. S1-S4. The convergence error is defined as where Fm denotes the secondary radiation force obtained using the m th mesh_size parameter. In all cases, the probe particle was placed 0.35 distance from the scatterer, this distance is approximately halfway between the node and antinode and therefore has a non-zero force. Figure S1. Convergence analysis results when the scattering particle is aligned with the pressure antinode and the probe particle is along the z direction. The convergence error is plotted as a function of the mesh size parameter (top graph) and as a function of the number of degrees of freedom (bottom graph). Figure S2. Convergence analysis results when the scattering particle is aligned with the pressure node and the probe particle is along the z direction. The convergence error is plotted as a function of the mesh size parameter (top graph) and as a function of the number of degrees of freedom (bottom graph). Figure S3. Convergence analysis results when the scattering particle is aligned with the pressure antinode and the probe particle is along the r direction. The convergence error is plotted as a function of the mesh size parameter (top graph) and as a function of the number of degrees of freedom (bottom graph). Figure S4. Convergence analysis results when the scattering particle is aligned with the pressure node and the probe particle is along the r direction. The convergence error is plotted as a function of the mesh size parameter (top graph) and as a function of the number of degrees of freedom (bottom graph).
Except for the last case ( Figure S4), the convergence is fast, and no significant difference between the uniform and non-uniform mesh can be observed, as far as the number of degrees of freedom is concerned. For non-uniform mesh, for all four cases the error is below 1%, when mesh_size > 35, meaning that the particles are meshed between 0.0143 < mesh < 0.0286 and the other domains between 0.0238 < mesh < 0.0476 .
For the uniform mesh, the error is below 1% for all cases when mesh_size > 23, meaning that for all domains the discretization size is 0.0217 < mesh < 0.0435 .
As for the 2D case, increasing the number of degrees of freedom in the same order of magnitude as for the 3D case would require an inefficiently large mesh, we decided to use a quartic discretization in this case compared to the cubic discretization of the 3D case. The quartic discretization allows expansion of a medium-sized mesh into a large number of degrees of freedom using fourth order approximation of the solutions over each mesh element. Moreover, as for the 3D case no significant difference was observed between a uniform and non-uniform mesh, for the 2D investigations we only applied a simple uniform mesh with the same characteristics as before.
These results are summarized in Figures S5 and S6. Figure S5. Convergence analysis results when the scattering particle is aligned with the pressure node and the probe particle is along either the r or z direction. The convergence error is plotted as a function of the mesh size parameter. Figure S6. Convergence analysis results when the scattering particle is aligned with the pressure antinode and the probe particle is along either the r or z direction. The convergence error is plotted as a function of the mesh size parameter.
In both cases an extremely fast convergence can be observed irrespective of the particle position: the error is already less than 0.01% when the mesh scale parameter is 20. For the 2D case, the relationship between the mesh size parameter and the number of degrees of freedom is summarized in Table S1. Although the methods are not applicable for simulations of touching spheres, we show that the separation distance (the surface-to-surface distance) of the two spheres can be arbitrary low. In Figure S7, secondary radiation force results for PS particle in water, where the scatterer is at the antinode can be seen, showing separation distances as low as 0.001λ with successful simulation. Figure S7. Secondary radiation force when the scattering particle is aligned with the pressure antinode and the probe particle is along the z direction. The distance in this case corresponds to surface-to-surface distance of the particles.