Critique on "Volume penalization for inhomogeneous Neumann boundary conditions modeling scalar flux in complicated geometry"

Numerical simulation of multiphysics problems within complex domains has garnered much interest in the past couple of decades. In the seminal work by Angot et al. [1], the authors describe a simple approach for simulating the incompressible flow over obstacles by applying an additional forcing term to the governing equations. In [1], this volume penalization (VP) methodology (also known as the Brinkman penalization method) was used to impose no-slip Dirichlet boundary conditions at the obstacle interface. Due to the simplicity of its formulation and implementation, the VP technique has been successfully applied to study a variety of fluid-structure interaction problems, including but not limited to water entry/exit [2], wave energy conversion [3, 4], aquatic locomotion [5, 6], fluttering instabilities [7], and flapping flight of insects [8, 9]. In all of these applications, the Dirichlet boundary condition formulation of the VP method was used. In the past few years penalization methods for Neumann and more general Robin boundary conditions have been proposed, although the analysis of such techniques is still an active area of research [10, 11, 12, 13, 14]. Kadoch et al. [10] extended the Dirichlet boundary condition VP formulation of Angot et al. [1] to allow for the imposition of homogeneous Neumann boundary conditions. Independently within the context of distributed Lagrange multipliers based fictitious domain method, Doostmohammadi et al. [15] informally described a way to enforce homogeneous flux boundary conditions on an interface by simply setting the thermal conductivity to zero within the obstacle. Sakurai et al. [14] recently developed a flux-based VP framework for imposing inhomogeneous, spatially constant Neumann boundary conditions on the boundary of a penalization region, which formally extended the methodology of Kadoch et al. [10]. This extension enables the simulation of more complex problems within the VP framework, such as flux-driven thermal convection in irregular domains. In the flux-based VP approach of Sakurai et al., the diffusion coefficient of the governing equation is modified and an additional forcing term is applied near the interface in order to impose the desired flux value on the boundary. This provides a simple and efficient way to impose flux boundary conditions on embedded interfaces. Through empirical testing of the penalized Poisson equation, Sakurai et al. [14] conclude that their method degrades to first-order accuracy if the embedded interface is not grid-aligned/grid-conforming despite the use of second-order finite differences. They also conclude that if two interfaces are considered, grid-aligned or otherwise, and a different flux boundary condition is imposed on each of them, then the method also degrades to first-order spatial accuracy. However, the method is second-order accurate for grid-aligned interfaces if the same (spatially constant) Neumann boundary condition values are considered. In this letter, we provide counter-examples to demonstrate that it is possible to retain second-order accuracy using Sakurai et al.’s method, even when different flux boundary conditions are imposed on multiple interfaces that do not conform to the Cartesian grid. We consider both continuous and discontinuous indicator functions in our test problems. Both indicator functions yield a similar convergence rate for the problems considered here. We also find that the order of accuracy results for some of the cases presented in