PaperlayerView exact revision on arXiv

Source-linked AI summary

An unstructured finite-volume Helmholtz method with perfectly matched layers for heterogeneous two-phase acoustics

Chuanchao Xu, Jun Liu, Jan-Helge Dörsam, Mario Kupnik, Tomislav Maric, Dieter Bothe

arXiv:2610.01593v1physics.flu-dyncs.NAmath.NAphysics.comp-ph

TL;DR

Time-harmonic acoustics in gas–liquid media must handle discontinuous material properties, and previous approaches had not combined VOF, cell-centred polyhedral finite volumes, and PML truncation. This paper integrates those elements and verifies the formulation across wave, transmission, radiation, and force studies, including aligned-refinement convergence for the layered case.

  • Problem

    Gas–liquid acoustics combines discontinuous material properties with a need for open-domain truncation, but prior methods had not integrated VOF, cell-centred polyhedral finite volumes, and PMLs.

    Source: §1 · pp. 2–3

  • Method

    The paper combines a VOF-based two-phase material representation, cell-centred unstructured finite volumes, and complex-coordinate PMLs in a real-valued block formulation.

    Source: §1 · p. 3

  • Results

    With aligned interface and PML-start refinement, the layered case’s whole-domain complex-pressure error decreases monotonically to 2.657 × 10^-4; the studies also report convergence for piston far-field pressure and rigid-sphere force verification.

    Source: §5 · p. 35, §6 · p. 36

  • Takeaways & Limitations

    The verification studies quantify the accuracy of the coupled interface–PML solution and assess radiation-pressure postprocessing across the tested cases.

    Source: §5 · p. 35

    The interface flux remains a one-field closure rather than an exact sharp cut-interface quadrature, so its accuracy is assessed under interface displacement and mesh refinement.

    Source: §3.2 · p. 12, §5 · p. 35

Abstract

from arXiv · show

This paper presents a cell-centred unstructured finite-volume method for time-harmonic two-phase acoustics. A volume-fraction representation supplies density and acoustic compressibility. Face fractions average available neighbouring planar interface cuts independently of velocity direction, retaining the established one-field face-density closure. Cartesian complex-coordinate stretching provides perfectly matched layers (PMLs). Real-imaginary splitting yields a real-valued block system, discretized with corrected non-orthogonal fluxes and consistent boundary contributions. Verification covers homogeneous waves, layered gas-liquid transmission with PML truncation, baffled-piston radiation with Kirchhoff far-field reconstruction, and resolved rigid-sphere radiation forces. Homogeneous-wave pressure converges at approximately second order on orthogonal meshes and on meshes with non-orthogonal interiors and orthogonal boundary cells. Reconstructed velocity converges at approximately second order on orthogonal meshes and with an observed order of about 1.7 on the latter mesh family. Under aligned refinement, the layered case's relative whole-domain complex-pressure $L_2$ error decreases to $2.657\times10^{-4}$. The resolved-sphere force differs from the Gorkov prediction by at most 1.5% over the tested Rayleigh size range. The results quantify the accuracy and current limitations of the heterogeneous Helmholtz-PML formulation.

1. Introduction

Time-harmonic acoustics in gas–liquid configurations must account for discontinuous material properties. This paper integrates a volume-fraction description, cell-centred unstructured finite volumes, and PML truncation in one frequency-domain method.

Motivation

  • Interface reflection and transmission affect the pressure and momentum-flux fields that determine an acoustically levitated drop’s equilibrium and stability.Time-harmonic formulations are used when the response at a prescribed frequency is sought, without resolving many acoustic periods.Source: §1 · p. 2

Research gap

  • Earlier finite-volume Helmholtz, interface-specific, and geometrical VOF methods addressed separate parts of the problem, but not their combination on cell-centred polyhedral meshes with PML truncation.The cited interface treatments use structured grids, while the cited geometrical VOF methods do not discretize the pressure Helmholtz formulation.Source: §1 · pp. 2–3

Scope

  • The model assumes linear, adiabatic acoustics about a quiescent base state with fixed interfaces, and excludes mean flow, viscosity, heat conduction, nonlinear acoustics, and interface motion over an acoustic period.A companion software manuscript addresses implementation and performance; this paper focuses on formulation and verification.Source: §1 · p. 3

Contribution

  • The method combines a VOF-based gas–liquid representation, cell-centred unstructured finite volumes, and PML truncation, with face fractions averaged from available adjacent-cell interface cuts.The shared face fraction is independent of velocity direction, and the complex-coordinate PML is expressed as a real-valued block system.Source: §1 · p. 3

2. Mathematical model

The mathematical model derives a pressure wave equation for spatially heterogeneous materials, then recasts it for time-harmonic fields and radiation-pressure evaluation. Complex-coordinate stretching supplies the PML formulation, whose real and imaginary pressure components form a coupled real-valued system.

2.1. Wave equation for heterogeneous media

  • Linearization about a quiescent heterogeneous base state with adiabatic closure yields a pressure wave equation whose density and sound-speed coefficients retain spatial dependence.The homogeneous wave equation follows when density and sound speed are taken as constants.Source: §2.1 · pp. 4–5

2.2. Helmholtz equation

  • The time-harmonic complex pressure satisfies a conservative Helmholtz equation with a = 1/ρ0 and β = 1/(ρ0c²); sharp-interface geometry is fixed and first-order capillary–acoustic coupling is neglected.The stated capillary-coupling assumption is appropriate away from capillary resonances when interfacial stress is small relative to acoustic stress.Source: §2.2 · pp. 5–6

2.3. Acoustic radiation pressure

  • For a homogeneous medium, time-averaged acoustic radiation pressure is evaluated from the first-order pressure and velocity fields after a second-order field expansion.The derivation assumes constant unperturbed density and sound speed.Source: §2.3 · pp. 6–7

2.4. Perfectly matched layer

  • Coordinate-wise complex stretching attenuates waves in homogeneous exterior PML regions; real–imaginary splitting yields coupled equations, while finite discretized layers retain truncation error.The physical-domain PML correction coefficients vanish, preserving the original heterogeneous operator there.Source: §2.4 · pp. 7–8

3. Unstructured FVM Discretization

The cell-centred polyhedral FVM combines a two-phase volume-fraction representation with face-based material coefficients and corrected fluxes. Its discretization carries interface geometry into the acoustic flux while accommodating non-orthogonal meshes, boundary conditions, and a real-valued block system.

3.1. Two-phase material representation

  • The one-field model linearly averages density and acoustic compressibility from the liquid volume fraction, and uses the same closure for face density.α = 1 in liquid, α = 0 in gas, with intermediate values in interface cells.Source: §3.1 · pp. 10–11
  • Shared face fractions combine available adjacent-cell PLIC cuts, with processor-face averages formed independently of volumetric-flux direction.The face fraction changes the material coefficient without donor selection or additional pressure unknowns.Source: §3.1 · pp. 10–11

3.2. Interface flux

  • The interface flux retains the one-field face coefficient and common face gradient, but does not impose sharp-interface conditions exactly in cut cells.Its accuracy therefore requires assessment under interface displacement and mesh refinement.Source: §3.2 · p. 12

3.4. Internal and boundary face gradients

  • Corrected face gradients split fluxes into implicit orthogonal and explicit non-orthogonal contributions, while boundary conditions supply fluxes and corresponding matrix or right-hand-side terms.The correction vanishes on orthogonal meshes; prescribed complex normal velocity is imposed through the Neumann contribution.Source: §3.4 · p. 13

3.5. Block system and correction sequence

  • Real-imaginary splitting produces a block system that is reassembled after explicit face-gradient corrections in each corrected-scheme sweep.Nominal second-order accuracy is expected for smooth coefficients and sufficiently regular meshes, while interface representation and reconstructed gradients can affect accuracy.Source: §3.5 · pp. 13–14

4. Verification

Verification spans homogeneous plane waves and layered gas–liquid transmission with directional PMLs. Refinement and damping studies quantify accuracy, including the sensitivity of reconstructed velocity to mesh deformation.

4.1. Computational workflow and AI assistance

  • The verification workflow used compiled solver results or reprocessed archived outputs, and the secondary-data package records plot and table reproduction inputs and scripts.The corresponding author directed the revision tasks and selected the included changes.Source: §4.1 · p. 14

4.2. Error measures

  • Observed order is estimated from successive mesh sizes; layered fields are probed at 1200 uniformly spaced points, while cell-field errors use volume weights.The layered sampling spans the complete domain, including endpoints and the PML.Source: §4.2 · pp. 14–15, §4.4 · p. 20

4.3. Homogeneous plane-wave baseline

  • The baseline tests oblique Dirichlet and aligned mixed Dirichlet–Neumann waves on orthogonal and smoothly deformed-interior meshes retaining an orthogonal boundary-cell layer.The interior deformation reaches maximum non-orthogonality of 32.8°–37.2°; two correction sweeps are used.Source: §4.3 · pp. 15–16, Figure 3 · p. 15, Figure 4 · p. 16
  • On orthogonal meshes, pressure order is 2.00–2.01 and velocity is approximately second order; deformed interiors yield pressure orders 2.04/1.94 and velocity orders 1.67/1.68.When boundary-adjacent cells are also deformed, finest pressure orders fall to 1.03 and 0.82, mixing interior-flux and boundary-condition errors.Source: Figure 5 · p. 16, Table 1 · p. 17, §4.3 · p. 18

4.4. Layered gas–liquid interface with PML

  • For the gas-to-liquid layered case, whole-domain EP is 1.197 × 10^-3, including the interface cell and PML attenuation.The comparison covers the reflected and transmitted fields together; figures show real and imaginary components and pressure-amplitude decay through the PML.Source: Figure 7 · p. 19, Figure 8 · p. 19, §4.4 · p. 20
  • For reverse water-to-air incidence, R = −0.999451 and whole-domain EP = 4.425 × 10^-3, with almost complete phase-reversed reflection.The transmitted pressure in air is much smaller, and numerical and analytical complex pressures agree closely over the domain.Source: §4.4 · p. 20, Figure 9 · p. 21

4.4.1. PML damping and mesh sensitivity

  • Increasing σmax from 10^4 to 5 × 10^5 s^-1 reduces EP from 5.188 × 10^-1 to 1.197 × 10^-3.At the selected damping, the real- and imaginary-component errors are 1.303 × 10^-3 and 1.122 × 10^-3, respectively.Source: §4.4.1 · p. 21, Figure 10 · p. 22, Table 2 · p. 22
  • EP falls monotonically from 1.997 × 10^-2 at Nx = 560 to 2.657 × 10^-4 at Nx = 8960 under aligned refinement.Both component errors also decrease monotonically over this sequence.Source: Figure 11 · p. 23, §4.4.1 · p. 23, Table 3 · p. 23

4.5. Controlled comparison of face-area definitions

The face-area comparison finds essentially unchanged pressure accuracy alongside reduced decomposition dependence for geometric averaging. The baffled-piston study assesses near- and far-field pressure against analytical references, with errors decreasing under mesh refinement and the angular pattern following the analytical result.

4.5. Controlled comparison of face-area definitions

  • Geometric averaging uses the same prepared mesh and initial interface while retaining the pressure flux, interpolation, quadrature, boundary conditions, PML coefficients, and velocity postprocessing.The comparison tests planar and spherical cases against an independent spherical-series pressure reference, with no PML in the sphere comparison.Source: §4.5 · pp. 23–24
  • The selected and original face-area methods agree to roundoff in the layered checks, with a largest relative field difference of 1.95 × 10−12.The checks cover 44 runs on one/four ranks and 32 runs on one/eight ranks.Source: §4.5 · p. 24
  • On the finest tested sphere mesh, geometric averaging changes EP from 2.9746% to 2.9753% and reduces D8 from 2.64 × 10−6 to 4.92 × 10−15.Both selections retain 3,041,280 matrix entries; the comparison supports reproducibility, not improved force accuracy or general curved-interface convergence.Source: Table 4 · p. 24, §4.5 · pp. 24–25
  • The comparison does not establish improved force accuracy or a general convergence order for curved interfaces.The reported evidence supports reproducibility with essentially unchanged pressure accuracy for the tested cases.Source: §4.5 · p. 25

4.6. Baffled-piston radiation

  • The baffled-piston case uses a rectangular PML around the outer radial and upper boundaries, with the piston edge fixed during refinement.The numerical near- and far-field pressures are compared with analytical references.Source: §4.6 · p. 25, Figure 12 · p. 25, Figure 13 · p. 27
  • The far field is reconstructed from numerical pressure and its gradient on a closed Kirchhoff–Helmholtz surface enclosing the piston source.The upper contour is reflected across the rigid baffle and revolved for the three-dimensional quadrature.Source: §4.6 · pp. 26–27
  • Near-field E2 decreases from 2.59 × 10−2 at 20 cells per wavelength to 2.58 × 10−3 at 80 cells per wavelength.The piston edge remains fixed across meshes, and the reconstruction surface lies in the undamped physical domain, meeting the PML onset at its endpoints.Source: §4.6 · pp. 25, 29, Table 5 · p. 27
  • Far-field pressure-magnitude error decreases from 2.75 × 10−2 at 20 cells per wavelength to 8.04 × 10−3 at 80 cells per wavelength.The reported angular pressure pattern follows the analytical solution.Source: Figure 14 · p. 28, Table 6 · p. 29, §4.6 · p. 29

4.7. Rigid-sphere radiation force in a standing wave

The rigid-sphere study evaluates radiation-stress force against the Gorkov prediction in a closed standing-wave cavity. The tested radius range recovers the expected Rayleigh scaling, and the position sweep reproduces the force lobes, signs, and zero crossings.

4.7. Rigid-sphere radiation force in a standing wave

  • For the Rayleigh-regime reference, the rigid-sphere limit uses f1 = f2 = 1 and Φ = 5/6.The Gorkov potential provides the point-particle comparison for the resolved-sphere force.Source: §4.7 · p. 30
  • The force is computed from time-averaged radiation-stress traction on the resolved sphere, with the normal pointing from fluid into sphere.The momentum-flux term is retained in the discrete traction expression although its contribution is zero at an exactly sound-hard wall.Source: §4.7 · p. 30
  • The reference case is a closed, rigid-walled cylindrical cavity with no active PML; only transducer excitation is prescribed, while pressure and scattered field around the sphere are computed.The cavity has height 7.64 mm, and the sphere-adjacent mesh resolution is parameterized by hs/a.Source: §4.7 · p. 30, Figure 16 · p. 31
  • Refining from hs/a = 1/8 to 1/16 changes force by 2.665%, while refinement from 1/16 to 1/32 changes it by 0.199%.The medium and fine results differ from the Gorkov value by 0.104% and 0.095%, respectively; the medium level is used for parameter sweeps.Source: Table 7 · p. 31, §4.7 · p. 32
  • The maximum difference from the Gorkov prediction is 1.485% over radii 30–150 µm, corresponding to 0.0139 ≤ ka ≤ 0.0694.The numerical force follows the expected a^3 scaling within the Rayleigh regime.Source: §4.7 · p. 32
  • For a 60 µm sphere, the position sweep reproduces both force lobes, their signs, and two analytical zero crossings; the RMS force difference is 1.288% of peak Gorkov force.The maximum peak-normalized difference is 3.415% at the first analytical zero, where the numerical residual is 0.0633 fN.Source: §4.7 · p. 32, Figure 19 · p. 34

5. Discussion

The discussion finds convergence and radiation benchmarks support the method within tested conditions, while identifying limits in coefficient treatment and mesh geometry. PML accuracy depends on sufficient attenuation, and the material closure does not establish arbitrary-interface accuracy.

5. Discussion

  • The velocity-independent face-fraction closure yields essentially unchanged tested pressure errors and reduces decomposition dependence to roundoff, but does not impose exact cut-cell interface conditions.It is a one-field closure; second-order accuracy for arbitrary subcell interface positions is not claimed, and sharp-interface flux treatments are outside scope.Source: §5 · p. 35
  • Local PML coefficients depend on the damping tensor and acoustic reaction coefficient, preserving the real-valued equations and block structure in both physical and absorbing regions.The layered damping study indicates that outgoing waves must be sufficiently attenuated before reaching the outer boundary.Source: §5 · p. 35
  • Pressure converges at approximately second order on orthogonal and non-orthogonal-interior meshes, but boundary-skewed diagnostics lose pressure order and arbitrary polyhedral geometry is not covered by a general second-order claim.Reconstructed velocity converges more slowly on the non-orthogonal-interior mesh family.Source: §5 · p. 35
  • Aligned refinement monotonically reduces the layered case's whole-domain complex-pressure error, while piston tests show far-field pressure convergence and the angular radiation pattern.The layered damping study also shows reduced accuracy when outgoing waves are not sufficiently attenuated before the outer boundary.Source: §5 · p. 35
  • The rigid-sphere study recovers Rayleigh a^3 scaling and force phase and sign changes, with sub-percent medium-to-fine change and a small sound-hard-surface momentum-flux contribution.The integrated force is controlled by the radiation-pressure field rather than an unresolved normal-velocity term.Source: §5 · p. 35

6. Conclusions

The conclusions present an unstructured Helmholtz–PML finite-volume method for heterogeneous two-phase acoustics and report benchmark accuracy for waves, transmission, and radiation. They identify additional acoustic-loss and capillary-coupling physics as possible extensions.

6. Conclusions

  • The method combines volume-fraction material representation, a real-valued split PML block system, and geometric face averaging in an unstructured Helmholtz finite-volume formulation.Pressure converges at second order on orthogonal and non-orthogonal-interior meshes; velocity converges at about order 1.7 on the latter family.Source: §6 · p. 36
  • With the interface and PML start aligned during refinement, the layered-case relative whole-domain complex-pressure error decreases monotonically to 2.657 × 10^-4 at Nx = 8960.The sequence demonstrates convergence of the coupled interface–PML solution.Source: §6 · p. 36
  • 0.258% near-field error and 0.804% far-field pressure-magnitude error validate piston radiation; sphere force differs from Gorkov by at most 1.485% over tested Rayleigh sizes.The sphere position-sweep peak-normalized RMS difference is 1.288%.Source: §6 · p. 36
  • Future work could add bulk viscosity and solid-wall dissipation, then couple the first-order acoustic field to streaming and interface-displacement models.A further extension could include first-order capillary–acoustic interaction through a linearized Young–Laplace condition.Source: §6 · pp. 36–37

Code and data availability

The implementation and secondary verification data are maintained in a repository, with provenance records mapping figures and tables to sources. Raw simulation fields and meshes are not included in the secondary-data package, and public release of the archive remains pending.

Code and data availability

  • The implementation and case workflows are maintained in the TwoPhaseFlow repository on the feature/acousticLevitation development branch.The submission cases specify the plicAverage area-fraction selection and legacy flux entry.Source: Code and data availability · p. 37
  • The implementation is identified at commit bb09dbfb34c1, while the deposited archive's DOI registration and public release remain pending.Version 1.0 is deposited in TUdatalib, with files identified by the checksum manifest and deposit record.Source: Code and data availability · pp. 37–38
  • Figure and table inputs are mapped in PROVENANCE.md, and the secondary-data package contains sampled profiles, convergence metrics, integrated forces, references, dictionaries, and regeneration scripts.Raw simulation fields and meshes are not included; the package has a checksum manifest.Source: Code and data availability · pp. 37–38

Funding

The authors acknowledge DFG funding for their tandem project and state that the funder had no role in the study or publication decision.

Funding

  • The tandem project was funded by the DFG, project number 542327521; the funder had no role in study design, data work, manuscript preparation, or submission.Mario Kupnik and Dieter Bothe acknowledge the support.Source: Funding · p. 38

CRediT authorship contribution statement

The statement assigns each author's contributions across the work. Xu and Maric have the broadest listed roles, spanning methodology, software, and writing-related work.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

During revision, OpenAI Codex assisted with manuscript and software-related work under the corresponding author's direction. The authors retain responsibility for the scientific content, interpretation, and final manuscript.

Loading 2610.01593v1…