Handling Neumann and Robin boundary conditions in a fictitious domain volume penalization framework

Sakurai et al. (J Comput Phys, 2019) presented a flux-based volume penalization (VP) approach for imposing inhomogeneous Neumann boundary conditions on embedded interfaces. The flux-based VP method modifies the diffusion coefficient of the original elliptic (Poisson) equation and uses a flux-forcing function as a source term in the equation to impose the Neumann boundary conditions. As such, the flux-based VP method can be easily incorporated into existing fictitious domain codes. Sakurai et al. relied on an analytical construction of flux-forcing functions, which limits the practicality of the approach. Because of the analytical approach taken in the prior work, only (spatially) constant flux values on simple interfaces were considered. In this paper, we present a numerical technique for constructing flux-forcing functions for arbitrarily complex boundaries. The imposed flux values are also allowed to vary spatially in our approach. Furthermore, the flux-based VP method is extended to include (spatially varying) Robin boundary conditions, which makes the flux-based VP method even more general. We consider several two- and three-dimensional test examples to access the spatial accuracy of the numerical solutions. The method is also used to simulate flux-driven thermal convection in a concentric annular domain. We formally derive the flux-based volume penalized Poisson equation satisfying Neumann/Robin boundary condition in strong form; such a derivation was not presented in Sakurai et al., where the equation first appeared for the Neumann problem. The derivation reveals that the flux-based VP approach relies on a surface delta function to impose inhomogeneous Neumann/Robin boundary conditions. However, explicit construction of the delta function is not necessary for the flux-based VP method, which makes it different from other diffuse domain equations presented in the literature.


Introduction
Partial differential equations (PDEs) in complex domains describe many natural and engineering processes. Examples include heat and mass transfer across melting/solidifying fronts, aquatic locomotion, cellular phenomena like cellular blebbing and cell crawling, flow in internal combustion engines or left ventricular assist devices, energy harvesting using wind turbines and wave energy converters, etc. In order to obtain meaningful solutions to PDEs, appropriate boundary conditions are required on the domain boundaries. Traditionally, body-fitted grid approaches, in which a complex domain is triangulated using sophisticated grid generation software, have been employed to solve PDEs numerically. Although body-fitted grid approaches allow imposing various types of boundary conditions accurately, they pose a serious challenge when the solution domain changes its topology over time. Issues like constant remeshing of the computational domain, the high aspect ratio of the elements, etc., limit the feasibility of body-fitted grid methods for modeling challenging moving domain problems.
indicator functions in the test problems considered in Sec. 5. It is observed that the continuous indicator function performs better (in terms of order of accuracy and uniformity of convergence rate) for imposing the spatially constant Neumann/Robin boundary condition, whereas the discontinuous one performs better for the spatially varying Neumann/Robin problem.
We also provide a formal derivation of the flux-based VP Poisson equation, which was not provided in Sakurai et al. [14], where the equation first appeared for the Neumann problem. The derivation reveals that the flux-based volume penalization method also uses a surface delta function to impose inhomogeneous Neumann/Robin boundary conditions. Interestingly, explicit construction of the delta function is not required in the flux-based approach, which is in contrast to the volume penalization approach of Ramière et al. [12]. We remark that on a formulation level the volume penalization approaches of Ramière et al., Kadoch et al., and Sakurai et al. (and the present work) are equivalent; minor differences in these works arise from the definition of the surface delta function. This insight is gained from Li et al. [18] who derived phase field-based diffuse domain equations satisfying Dirichlet, Neumann, and Robin boundary conditions. Li et al. used the method of matched asymptotic expansions to provide different diffuse domain approximations for the Neumann/Robin problem 1 ; these approximations differ in the way how surface delta function is defined.
Characteristic-based approaches to impose Neumann and Robin boundary conditions for the volume penalized PDEs have also been proposed in the literature; see, for example, Brown-Dymkoski et al. [19] and Hardy et al. [20] who used characteristic-based VP approach to model the energy transport equation satisfying Neumann and Robin boundary conditions in the context of compressible flows and low Mach formulation of compressible flows, respectively. The main limitation of the characteristic-based VP method is that it relies on having a time-derivative term in the PDE and as such cannot be applied to steadystate (i.e., having no temporal derivative term) elliptic equations. In addition to the volume penalization methods [12,13,14,19,20,21,22,23], other fictitious domain techniques have also been proposed to impose flux boundary conditions on embedded interfaces. Notable ones include the flux-correction technique (FCT) of Ren et al. [24], Wang et al. [25], and Guo et al. [26] and the direct forcing method of Lou et al. [27]. FCT is a predictor-corrector scheme and is implemented using the Lagrangian-Eulerian machinery of the IB method. In the prediction step of FCT, an intermediate scalar field is computed on the Eulerian grid, which in general does not satisfy the flux boundary condition on the interface defined by the Lagrangian markers. Next, in the correction step, a Lagrangian forcing term is computed either implicitly [25,26] or explicitly [24] that corrects the intermediate scalar field to satisfy the Neumann boundary condition. In an essence, FCT is a time-splitting approach (similar to the characteristic-based VP approach), which requires having a time-derivative term in the scalar transport equation. Therefore, unlike the flux-based VP method, FCT cannot be used for time-independent elliptic equations. In the direct forcing method, the scalar field near the interface is reconstructed locally using second-or third-degree polynomials in order to satisfy the flux boundary condition. This is achieved by identifying "forcing" points on the fictitious (solid) side of the interface, on which the reconstructed scalar field value is directly imposed. Direct forcing methods are also typically implemented as a predictor-corrector scheme, which avoids modifying the system of linear equations.
In the following sections, we first describe the continuous form of the volume penalized equations and thereafter describe the numerical construction of the flux-forcing functions. Finally, various test cases are considered in two-and three-spatial dimensions to access the accuracy of the numerical solutions.

The Neumann problem
Consider an irregular fluid domain Ω f embedded into a larger, regular computational domain Ω, as shown in Fig. 1. Define Ω \ Ω f = Ω s as the fictitious solid domain and n as a unit outward normal of the fluid-solid interface ∂Ω s (or ∂Ω f ). With q as the scalar quantity of interest, κ as the diffusion coefficient, and f as a satisfying inhomogeneous Neumann/flux boundary conditions on ∂Ω s − κ n · ∇q = g, to the entire computational domain Ω using the flux-based VP approach. The extended domain Poisson equation satisfying the inhomogeneous flux boundary conditions on the interface reads as Here, η is the penalization parameter, χ(x) is an indicator function whose value is 1 in the solid region and 0 in the fluid region, and f b = ∇ · (χβ) − χ∇ · β is an additional forcing term required to impose the flux boundary conditions on ∂Ω s . The vector-valued flux-forcing function β(x) is selected such that β ·n = −g on the interface. In the limit of η → 0, the solution to the volume penalized (VP) Poisson equation converges to the solution of non-penalized Poisson equation (Eqs. (1) and (2)). A formal derivation of Eq. (3) is provided in Appendix A. As noted in Thirumalaisamy et al. [15], the flux-based VP approach allows κ and g to vary spatially as well.

The Robin problem
Next, we consider the inhomogeneous Robin boundary conditions of the type on the fluid-solid interface ∂Ω s . Appendix B derives the flux-based VP Poisson equation for the Robin problem, which reads as In the equation above, the flux-forcing function satisfies the requirement of β · n = −g. The unit normal vector n appearing in the first term of Eq. (5) can be computed numerically using a signed distance function as explained later in Sec. 2.4. In our formulation, ζ, κ, and g are allowed to vary spatially.

Multiple interfaces and coupled volume penalized equations
The VP Poisson equations (Eqs. (3) and (5)) can also be generalized to handle multiple interfaces within the computational domain Ω. For some of these interfaces, Dirichlet boundary conditions may also be prescribed. Following Thirumalaisamy et al. [15], the generalized form of the VP Poisson equation satisfying Neumann and Dirichlet boundary conditions reads as (6) For the above equation to hold true, the computational domain Ω is assumed to consist of disjoint volumetric regions Ω d i (for i = 1, 2, . . . , D) and Ω n j (for j = 1, 2, . . . , N ), with imposed Dirichlet (q = q d i ) and Neumann (−κ ∇q · n j = g n j ) boundary conditions, respectively. Furthermore, the union of Ω d i and Ω n j regions defines the total solid domain, i.e., Ω s = Ω d 1 ∪ Ω d 2 ∪ · · · Ω d D ∪ Ω n 1 ∪ Ω n 2 ∪ · · · Ω n N . In Eq. (6) the indicator function χ n (x)(respectively, χ d (x)) is 1 if x ∈ Ω n (respectively, Ω d ) and 0 if x ∈ Ω \ Ω n (respectively, Ω d ). Note that Robin boundary conditions can be easily included in Eq. (6), as their form is very similar to the Neumann problem. We omit Robin boundary conditions in the generalized equation written above for brevity.
The volume penalization approach can also be extended to other governing equations that describe conservation of momentum, energy, species, etc. For example, the VP incompressible Navier-Stokes equations coupled to the flux-based VP advection-diffusion equation satisfying Neumann boundary condition reads as ∇ · u = 0, In the equations above, u(x, t) is the fluid velocity, u b (x, t) is the structure velocity, p(x, t) is the hydrodynamic pressure, f denotes the momentum body force, ρ(x) is the mass density, and µ(x) is the dynamic viscosity. The equation set (7)-(9) is written considering only a single interface in the domain; generalization to handle multiple interfaces is also possible following Eq. 6. We remark that in the context of fluid-structure interaction (FSI) problems, only velocity matching condition on the fluid-structure interface is required, i.e., u = u b on ∂Ω s is essential, whereas u = u b in Ω s is optional. In the volume penalization approach to FSI, both these conditions are imposed through the penalization term χ η (u b − u). Therefore, in the VP momentum equation (7), only Dirichlet boundary conditions have been considered.

Interface capturing
We use a signed distance function φ(x) to implicitly define the fluid-solid interface ∂Ω s . The scalar field φ(x) is defined to satisfy the following property: Moreover, the negative gradient of the signed distance function φ(x) gives the unit outward normal vector of the interface, i.e., n = −∇φ. The signed distance function can also be used to define the indicator function χ(x). In this work we use φ(x) to define two types of indicator functions: one is smooth and continuous and written as and the other one is discontinuous, which reads as In the Eq. (10) above, n smear ∈ R is the number of grid cells over which the indicator function is smoothed on either side of the interface and h is the grid cell size.

Discrete equations
We use second-order finite difference stencils to discretize the spatial derivative terms of cell-centered Poisson and face-centered momentum equations (Eqs. (3) and (7), respectively) on a Cartesian grid. Fig. 2 shows a schematic representation of a two-dimensional Cartesian grid cell, in which the velocity components are stored on edge centers (face centers in three spatial dimensions), whereas the transported variable q, the fluid pressure p, and the signed distance function φ are stored at the cell center. The computational domain Ω is discretized into N x × N y Cartesian grid cells with mesh spacing ∆x and ∆y in the x-and y-direction, respectively. In this work we use equal mesh spacing in the two directions, i.e., ∆x = ∆y = h. In what follows next, we primarily focus on the discretization of the VP Poisson Eq. (3) for the Neumann problem; details on the spatiotemporal discretization of the VP incompressible Navier-Stokes equations can be found in our prior works [28,6,4].
Referring to Fig. 2, let (i, j) denote the cell index, (i − 1 2 , j) denote the lower x edge index and (i, j − 1 2 ) denote the lower y edge index. Then the discretized form of the VP Poisson Eq. 3 in two spatial dimensions reads as in which and the right hand side term S i,j is given by Analogous discretization formulas can be written for the three-dimensional VP Poisson equation. In the discretized Eq. (12) written above, the indicator function χ and the diffusion coefficient κ are required at the edge centers; these properties are first defined at the cell centers and then interpolated onto the edge centers using a second-order accurate linear interpolation scheme. The flux-forcing function β(x) is also required at the edge centers; methods to construct β are discussed next.

Construction of flux-forcing functions
The vector-valued flux-forcing function β(x) plays a crucial role in imposing the desired inhomogeneous Neumann and Robin boundary conditions on the interface. In this section, we introduce three approaches to construct β(x), namely Approach A, B, and C. The three approaches are in increasing order of generality. While Approach A is specialized for the Neumann problem, Approaches B and C are equally applicable to the Robin problem.
3.1.1. Approach A: Analytical construction of spatially varying g Consider for a moment that the solution to the non-penalized Poisson equation with inhomogeneous Neumann boundary conditions on ∂Ω s is known. Denote the exact solution by q exact . If the flux-forcing function is taken to be of the form β(x) = κ ∇q exact (x), then it satisfies the requirement of β · n = −g on ∂Ω s . In practice the solution to the Poisson Eq. (1) with boundary condition (2) is sought and not known a priori. However, if an analytical approximation q to the exact solution q exact exists, such that n · ∇ q = n · ∇q exact on ∂Ω s , then can be prescribed as a flux-forcing function. Away from the interface, the approximation q can be close to or very different from q exact , depending upon whether a continuous or a discontinuous indicator function χ is used. We denote the analytical construction of β as Approach A. In component form, Approach A is written as Approach A was employed in Sakurai et al. [14] and Thirumalaisamy et al. [15] to demonstrate the feasibility of flux-based volume penalization method to solve PDEs with flux boundary conditions in complex domains, but is quite restrictive in practice as discussed next.

Approach B: Numerical construction of spatially constant g
Although Approach A allows for imposing spatially varying g values on the interface, approximating an analytical solution to the exact solution near the interface is a non-trivial task, especially if the interface is geometrically complex. However, if g is spatially constant, then constructing β is easy. This is achieved by taking β(x) = −g n(x), as it satisfies the requirement of β · n = −g on ∂Ω s . Now, recalling from Sec. 2.4 that the negative gradient of the signed distance function φ(x) is the continuous normal vector field n(x), the flux-forcing function can be constructed numerically for an irregular boundary as We denote the numerical construction of spatially constant g value on the interface as Approach B, which in component form is written as  3.1.3. Approach C: Numerical construction of spatially varying g As a generalization of Approach B, the flux-forcing function can be taken as However, Eq. (19) poses a challenge of extending the codimension-1 boundary condition function g defined over the interface to a codimension-0 function g(x) defined in the neighborhood of the interface. Although there are several ways to achieve this function extension (in absence of a constraint), in this work we follow a simple strategy of propagating the interfacial g values to the neighboring grid cells along the interface normal. More specifically, consider a fluid-solid interface ∂Ω s embedded into a Cartesian grid as shown in Fig. 3(A). The signed distance function φ(x) can be used to identify the grid cells through which the interface passes. Denote these grid cells as interface cells. Fig. 3(A) highlights two such interface cells: one whose cell center x lies in the fluid region and the other whose cell center x lies in the solid region. The normal vector of the interface cells is also known from the signed distance function: n = (−∇φ) and n = (−∇φ) . Next, the g value at the cell center of an interface cell is set equal to the closest interfacial g value: in which x ∂Ωs = x + φ n and x ∂Ωs = x + φ n are the closest points on the interface to the cell centers x and x , respectively. Note that the g function on the interface is prescribed and therefore, g(x ∂Ωs ) and g(x ∂Ωs ) are known a priori. In the next part of the algorithm, g and g values are propagated to the grid cells that are within a distance of n prop h to the interface cells along ± n and ± n directions, respectively. This procedure is pictorially described in Fig. 3(B). The number of grid cells n prop to which g values are propagated depends upon the choice of the indicator function χ-we will explore the effect of n prop on the solution accuracy in Sec. 5. Note that propagating g values along the normal directions may lead to a situation of conflict at a grid cell where two or more interface cell normals intersect. This situation is shown for the shaded cell in Fig. 3(B) where the two normals n and n intersect. For such cells, a g value with the larger modulus is chosen: We also considered an average and the minimum modulus of g at the conflicted cells; these choices however reduced the order of accuracy of the solution for the continuous/smoothed indicator function. Note that for the discontinuous indicator function, the g value at the conflicted cells does not matter much for the solution accuracy. This is because such cells are generally located far away from the interface where the discontinuous indicator function is already zero. Nevertheless, we always make use of Eq. (21) even for the discontinuous indicator function in this work. With g values defined at the cell centers, the component form of β reads as The propagation strategy of Approach C can also be implemented by solving a hyperbolic equation of the form The equation above can be integrated over a pseudo-time interval ∆τ that is directly related to the propagation distance. However, the test examples of Sec. 5 show that the method of g propagation described in Approach C is quite effective in imposing the spatially varying flux boundary conditions. Moreover, it does not require solving any additional partial differential equation.
Note that there can be other ways of extending the flux-forcing function in the vicinity of the interface as discussed at the beginning of this section. One straightforward approach is to extend the β function defined over the interface Ω s to a flux-forcing function valid near the interface using the top hat or a Gaussian belllike function, which we refer to as Approach D. However, as demonstrated in Appendix D, this particular function continuation approach does not produce satisfactory results; the numerical and actual solutions differ significantly and the numerical scheme does not converge under grid refinement in any norm. In contrast, Approach C produces the correct solution and a convergent numerical scheme. This also highlights the non-triviality in allowing spatially varying Neumann/Robin boundary conditions in the flux-based VP method.

Software
The flux-based volume penalization algorithms described here are implemented within the IBAMR library [29], which is an open-source C++ simulation software focused on immersed boundary and volume penalization methods with adaptive mesh refinement. The code and test cases presented in Sec. 5 are publicly available at https://github.com/IBAMR/IBAMR. IBAMR relies on SAMRAI [30,31] for Cartesian grid management and the AMR framework. Linear and nonlinear solver support in IBAMR is provided by the PETSc library [32,33,34]. All of the example cases in the present work made use of distributed-memory parallelism using the Message Passing Interface (MPI) library.

Results and discussion
In this section we discretely solve the volume penalized Poisson Eqs. 3 and 5 satisfying inhomogeneous Neumann and Robin boundary conditions, respectively, to assess the accuracy of the numerical solutions. We use the flexible GMRES (FGMRES) iterative solver with a tight relative residual tolerance of 10 −12 to solve the system of linear equations. The order of accuracy results presented here are computed only in the fluid domain and are determined based on the L 1 and L ∞ norm of the error (denoted E 1 and E ∞ , respectively) between the numerical and analytical 2 solutions. Since the VP method is expected to produce a non-uniform convergence rate under grid refinement because of the delta function formulation (see Appendices A and B for derivation), we curve-fit the error data and report the slope/convergence rate, denoted m and the coefficient of determination, denoted R 2 , in each case. Appendix F tabulates the error data. The spatial convergence rate of the error is shown for both continuous (denoted E ∞ c and E 1 c ) and discontinuous (denoted E ∞ d and E 1 d ) indicator functions. We consider two-and three-dimensional examples involving constant and spatially varying flux boundary conditions on ∂Ω s . In the test examples, the fluid region Ω f is embedded into a larger computational domain Ω with Dirichlet boundary conditions imposed on the external boundary ∂Ω of the domain. The computational domain is discretized into N × N and N × N × N grid cells for the two-and three-dimensional examples, respectively. The penalization parameter η is taken to be 10 −8 (Appendix E considers the effect of η on the convergence rate) and the diffusion coefficient κ is taken to be 1 for all of the tests. While imposing the Robin boundary conditions we take ζ to be 1 in the test examples. The numerical solutions are presented for Approach C and where applicable, results obtained from Approach C are compared against Approach A or B. Since Approach A constructs the flux-forcing function from the known solution to the problem, for a given indicator function χ, Approach A is expected to perform better than or at least as well as Approach B and C. This expectation is also confirmed from the tests that follow next.

Concentric circular annulus with spatially constant flux on the interface
We first consider the concentric circular annulus problem from Sakurai et al. [14] in which different inhomogeneous, but spatially constant, Neumann boundary conditions are specified on the two interfaces defining the annulus. The inner radius of the annulus is r i = π/4 and its outer radius is r o = 3π/4. The center of the annulus is positioned at (π, π). The circular annulus is embedded into a larger computational domain of extents Ω ∈ [0, 2π] 2 . The source term of the Poisson equation for this case is in which r = (x − π) 2 + (y − π) 2 and the Neumann boundary condition values on the two interfaces are taken to be dq dr r= π 4 = 3m and dq dr r= 3π The exact solution of this problem using the zero-mean condition Ω f rq(r) dr = 0 reads as The mean of the numerical solution in the fluid region is subtracted as a post-processing step to impose the zero-mean condition numerically. Relatively simple geometry and constant flux boundary conditions of this test problem allows for an analytical construction of β. Indeed, in [14,15], the flux-forcing function β was constructed analytically It is to be noted that d q dr reduces to 3m at r = π/4 and m at r = 3π/4, respectively. Since the flux boundary condition value is spatially constant on the interface, Approach B is also applicable for this test problem. The g value required for Approach B is m and −3m on the inner and outer interface, respectively. Figs. 4(A) and 4(B) show the order of accuracy of the solution as a function of mesh resolution for Approach A and B, respectively. As can be observed in Fig. 4(A), for Approach A, O(h 1.98 ) (respectively, O(h 1.94 )) convergence rate with an R 2 value of 0.99 (respectively, 0.96) in L ∞ (respectively, L 1 ) norm is achieved using the continuous indicator function. With the discontinuous χ, O(h 1.50 ) (respectively, O(h 1.78 )) convergence rate with an R 2 value of 0.98 (respectively, 0.99) in L ∞ (respectively, L 1 ) norm is  achieved. We note that the convergence rate using the continuous indicator function is better than the discontinuous function for Approach A. Looking at Fig. 4(B), it is seen that Approach B exhibits a very similar convergence rate as Approach A, but in contrast to Approach A, Approach B is more versatile as it requires only ∇φ(x) information, which can be constructed for any irregular boundary [17]. Figs. 4(C) and 4(D) compare the analytical and numerical solutions along x− and y−direction, respectively. Numerical solutions using Approach B and discontinuous indicator function are presented. As can be observed in the figures, an excellent agreement is obtained between the analytical and numerical solutions. Next, we solve this problem using Approach C. As can be seen in Fig. 5(A), Approach C also exhibits a very similar convergence rate as Approach A and B when the discontinuous indicator function is used, whereas the convergence rate is reduced when the continuous function is employed. Specifically, O(h 0.91 ) (respectively, O(h 1.22 )) convergence rate with an R 2 value of 0.79 (respectively, 0.82) in L ∞ (respectively, L 1 ) norm is obtained using the continuous χ. However, Approach C is the most general one, since it can be used for imposing spatially varying flux values as demonstrated in later examples.
The present example is also used to study the effect of the number of propagation cells n prop on the solution accuracy for the continuous indicator function. The results are shown in Fig. 5(B), in which it can be observed that 2 cells on either side of the interface are sufficient for propagating g values for a fixed number of n smear cells; the error norms are mostly affected by the n smear choice. Based on the results of this problem, we choose n smear = 1 and n prop = 2 for the continuous masking function, unless otherwise stated. For the discontinuous indicator function also, we use n prop = 2 for the remainder of the problems (although n prop = 1 is also sufficient).

Spatially varying flux values along complex interfaces
In this section, we assess the accuracy of the numerical solution for spatially varying flux values using a manufactured solution of the form q exact (x) = sin(x) sin(y).
Inhomogeneous Neumann boundary conditions g(x) = −κ n · ∇q exact are imposed on the fluid-solid interface ∂Ω s , whereas Dirichlet boundary conditions are imposed on the external boundaries of the computational domain, i.e., q| ∂Ω(x) = q exact (x). Note that g(x) varies spatially, and therefore, Approach B is not applicable for this test. Eq. (28)   Our prior work [15] considered Approach A for constructing the flux-forcing function for this problem, wherein q exact was used to define β = κ ∇q exact . Although not feasible in practice (solution is unknown), Approach A results in second-order convergence rate of the numerical solution for this problem; see Appendix C. Next, we solve the same problem using Approach C. Fig. 6 presents the numerical solution and its convergence rate as a function of grid resolution. For the hexagram case, at least O(h 0.56 ) accuracy is achieved using the continuous indicator function whereas at least O(h 0.78 ) accuracy is achieved using the discontinuous indicator function. Similarly, for the egg case, at least O(h 0.67 ) accuracy is achieved using the continuous χ and at least O(h 0.54 ) accuracy is obtained using the discontinuous indicator function. Lastly, for the x-cross geometry, Approach C exhibits at least O(h 0.67 ) accuracy using the continuous indicator function and at least O(h 0.85 ) accuracy with the discontinuous one. As noted in the previous section also, Approach C with the discontinuous indicator function is able to achieve a better convergence rate than with the discontinuous one. We remark that for Approach C, the reduction in accuracy (when compared to Approach A) is attributed to the codimension-0 extension of the spatially varying g function in the neighborhood of the interface. Nevertheless, Approach C is able to impose spatially varying flux values on a complex interface (sharp corners, etc.) and the solution accuracy is also reasonable. Later in Sec. 5.6 we demonstrate that smoothing of geometric features like sharp corners improves the accuracy of Approach C further.

Constant and spatially varying flux on three-dimensional interfaces
In this section, we consider two complex geometries in three spatial dimensions: a sphere and a torus. These geometries are embedded into a larger computational domain of extents Ω ∈ [0, 2π] 3 .
For the spherical geometry, the manufactured solution is taken to be The comparison between the numerical and analytical solutions, as well as the spatial convergence rate of E 1 and E ∞ error norms using Approach B and C are shown in Fig. 7. As can be observed in the figure, the numerical solution is in excellent agreement with the exact solution. Largely O(h 2 ) convergence rate is obtained for this example using Approach B with both continuous and discontinuous indicator functions. Approach C also yields the same order of accuracy with the discontinuous indicator function as Approach B. However, the convergence rate using the continuous indicator function is between 0 and 1 for Approach C. Clearly, the discontinuous indicator function performs better than the continuous one for Approach C. Second-order convergence rate is also obtained with Approach A, when β = κ ∇q exact is used for the spherical geometry (data not shown for brevity). Better performance of Approach A compared to Approach B and C is expected, as mentioned in the beginning of Sec. 5.
For the next three-dimensional test example, a solid torus is embedded in a computational domain of extents Ω ∈ [0, 2π] 3 as shown in Fig. 8(A) and the fluid region is taken outside of the torus. We consider a manufactured solution of the form q exact (x) = − cos(x) cos(y) cos(z).
(1) to generate the required source term f (x). On the toroidal interface, spatially varying inhomogeneous Neumann boundary conditions are imposed, whereas on the external domain boundary inhomogeneous Dirichlet boundary conditions using the exact solution are imposed. We solve this test problem using Approach C and the results are shown in Fig. 8. As can be observed in Fig. 8, at least O(h 1.25 ) is achieved with the discontinuous indicator function and at least O(h 0.63 ) is achieved with the continuous indicator function.

Spatially constant Robin boundary condition on two-dimensional interfaces
We consider the concentric circular annulus problem of Sec. 5.1 with the same exact solution q exact (r) and source term f (r), as written in Eqs. (26) and (24), respectively. Plugging the exact solution into the Robin boundary condition Eq. 4 yields a spatially constant g value for the inner and outer interface, respectively.
The VP Poisson Eq. 5 is solved using Approach B and C for this problem. The numerical solution compared against Eq. 26 using Approach C is shown in Fig. 9(A); an excellent agreement is obtained. We also present the convergence rate for Approach B and C in Fig. 9. As can be seen in Fig. 9, the convergence rates obtained by using the discontinuous indicator function for Approach B and C are quite close to what we had obtained in Sec. 5.1. For the continuous indicator function, we obtain approximately second-order accuracy with Approach B and at least O(h 0.71 ) accuracy with Approach C. As also observed in Sec. 5.1, the discontinuous indicator function performs better than the continuous one with Approach C; the reverse is true for Approach B.

Spatially varying Robin boundary condition on a complex two-dimensional interface
In this section, we assess the accuracy of Approach C for spatially varying Robin boundary conditions on a complex two-dimensional interface. A hexagram geometry is embedded into a computational domain of extents Ω ∈ [0, 2π] 2 , as considered in Sec. 5.2. The fluid is occupied between the computational domain boundary ∂Ω and the fluid-solid interface ∂Ω s . The same manufactured solution as written in Eq. 28 is considered here; this solution yields spatially varying g values when plugged into the Robin boundary condition Eq. (4). We solve the VP Poisson Eq. 5 using Approach C. The numerical solution and the spatial convergence rate of the error norms are presented in Fig. 10. As observed in the figure, at least O(h 0.71 ) accuracy is achieved with the continuous indicator function and at least O(h 0.72 ) accuracy is achieved with the discontinuous indicator function.

Effect of smoothing geometric features on the convergence rate of Approach C
In this section, we study the effect of sharp geometric features, such as corners on the convergence rate of Approach C. We consider a slightly modified version of the hexagram interface, which was considered earlier in Secs. 5.2 and 5.5. Instead of the sharp corners, in this case, we embed a hexagram having smooth exterior corners (see Fig. 11(A)) in a larger Cartesian domain of extents Ω ∈ [0, 2π] 2 and solve the Neumann/Robin problem using Approach C. The order of accuracy results for spatially varying Neumann and Robin boundary conditions using the discontinuous indicator function are presented in Fig. 11. For the sharp corner geometry case with spatially varying Neumann boundary conditions (as shown in Fig 6(A)), the convergence rates were O(h 0.78 ) and O(h 0.84 ) in L ∞ and L 1 norm, respectively. In contrast, with smooth corners, the convergence rates are O(h 0.95 ) and O(h 1.08 ) in L ∞ and L 1 norm, respectively. A similar trend is obtained when spatially varying Robin boundary conditions are considered: For the sharp geometry case (as shown in Fig. 10), the convergence rates were O(h 0.71 ) and O(h 1.00 ) in L ∞ and L 1 norm, respectively. With smooth corners, the convergence rates are O(h 1.00 ) and O(h 1.26 ) in L ∞ and L 1 norm, respectively. This test demonstrates that the convergence rate of Approach C also depends upon local geometric features.

Application to free convection problem
Finally, we consider steady natural convection in a concentric annulus. A constant heat flux Q is imposed on the inner cylinder of radius r i and a fixed temperature T o is maintained on the outer cylinder of radius r o . The concentric annulus is embedded into a larger computational domain of extents Ω ∈ [−2.56, 2.56] 2 , as shown in Fig. 12. This example was studied in Yoo [35] using a body-fitted grid approach and more recently, it has been used to validate the IB/FCT relying on time-splitting approach to handle the flux boundary conditions on embedded interfaces [24,25,26]; see Introduction Sec. 1 for a brief discussion on IB/FCT.  We solve the volume penalized advection-diffusion equation for the temperature field coupled to the volume penalized incompressible Navier-Stokes equations (Eqs. (7)-(9)) in a non-dimensional form. The nondimensional quantities are defined as: dimensionless temperature Φ * = k(T −T 0 )/(QL), velocity u * = uL/α, time t * = tα/L 2 , and position x * = x/L. Here, k is the thermal conductivity, α is the thermal diffusivity α = k/(ρc p ), ρ is the density, c p is the heat capacity at constant pressure, and L = r o − r i is the annulus thickness. For this case we consider r * i = r i /L = 1 and r * o = r o /L = 2. Dropping the * superscript from the non-dimensional quantities, the system of non-dimensional equations reads as Here, e y = (0, 1) is a unit vector in the y−direction and the continuous indicator functions χ d and χ n are defined to be In the above, φ d (x) and φ n (x) are the signed distance functions for the Dirichlet (outer cylinder) and Neumann (inner cylinder) boundary, respectively. The Rayleigh number Ra = GγQL 4 /(kαν) and the Prandtl number Pr = ν/α are the two main non-dimensional parameters that characterize buoyancy-driven flows; these parameters are seen in the right-hand side of the non-dimensional momentum Eq. (31). Here, γ is the coefficient of thermal expansion, G is the gravitational constant, and ν is the kinematic viscosity. The flux-forcing function β imposing the constant flux boundary condition on the surface of the inner cylinder, n i · ∇Φ = 1, is constructed using Approach B and C for this problem. Here, n i is the unit outward normal vector of the inner cylinder. On the outer cylinder homogeneous Dirichlet boundary condition, Φ = 0, is imposed through the last term of Eq. (33). Periodic boundary conditions are used on the external domain boundaries. Two Rayleigh numbers Ra = 5700 and Ra = 5 × 10 4 are considered for this problem. The same Prandtl number Pr = 0.71 is used for the two cases. The computational domain Ω is discretized by a uniform Cartesian grid of size 256 × 256. The penalization parameters η d and η n are taken to be 10 −8 . We treat the convective and the advective terms of Eqs. (31) and (33) explicitly, whereas the rest of the terms are treated implicitly. The implicit treatment of volume penalization terms in Eqs. (31) and (33) allows us to use a relatively large time step sizes of ∆t = 10 −4 and ∆t = 5 × 10 −5 for Ra = 5700 and Ra = 5 × 10 4 cases, respectively. In contrast, Sakurai et al. [14] used a time step size of ∆t = 10 −6 for these two cases as they employed an explicit Euler time marching scheme. More details on the second-order accurate spatial discretization and time-stepping scheme employed in the fluid solver can be found in our prior works [28,36].
To compare our results with those reported in [14] we plot the steady-state temperature distribution on the left half of the inner cylinder for both Ra cases in Fig. 12. In the figure, the polar angle Θ = 0 • starts from the top position (x, y) = (0, 1) of the inner cylinder and ends at its bottom position (x, y) = (0, −1), where Θ = 180 • . As observed in the figure, the numerical results obtained using both Approach B and C are in excellent agreement with those reported in [14] who used Approach A for constructing β. Sakurai et al. compared their numerical results with Yoo [35] and Ren et al. [24]; comparison with Yoo and Ren et al. is therefore omitted in Fig. 12 in the interest of clarity. We also present the steady-state temperature field in the whole annular domain at the two Rayleigh numbers in Fig. 12.

Conclusions
In this work, we proposed a numerical technique for constructing flux-forcing functions for the flux-based VP method introduced by Sakurai et al. We also extended the flux-based VP approach to include Robin boundary conditions. Our method of flux-forcing functions is more general than the analytical approach (denoted Approach A in this work) of Sakurai et al. and requires only a signed distance function to construct the flux-forcing function. Two numerical-based approaches were presented for constructing flux-forcing functions: Approach B for imposing spatially constant and Approach C for imposing spatially varying (as well spatially constant) Neumann/Robin boundary conditions. Within Approach C we extended the (spatially varying) codimension-1 g function to the neighborhood of the interface using a simple propagation strategy. We considered several two-and three-dimensional Poisson problems in complex domains to assess the accuracy of the numerical solutions. Results were presented for both continuous and discontinuous indicator functions. For Approach B, largely O(h 2 ) accuracy is observed using the continuous indicator function. Between O(h 1 ) and O(h 2 ) convergence rate is observed for Approach B with the discontinuous indicator function and a similar convergence rate is observed for Approach C with the discontinuous indicator function when it is used for solving the constant Neumann/Robin boundary condition problem. For spatially varying boundary conditions, Approach C using the discontinuous indicator function exhibits close to O(h 1 ) convergence rate; the accuracy of the method is further improved by smoothing the sharp geometric features. However, Approach C using the continuous indicator function exhibits a convergence rate between O(h 0 ) and O(h 1 ). Based on our results, we recommend using the discontinuous indicator function with Approach C to impose spatially varying Neumann/Robin boundary conditions and using the continuous indicator function with Approach B to impose spatially constant Neumann/Robin boundary conditions. Finally, we used Approach B and C to study the flux-driven thermal convection problem in a concentric annulus and compared our results against the literature. An excellent agreement was obtained. We also provided formal derivation of the flux-based volume penalized Poisson equations in strong form for both Neumann and Robin problems. The formulation shows that an explicit construction of the delta function is not necessary for the flux-based VP method, which makes it different from other diffuse domain equations presented in the literature. In this section, we derive Eq. (3) by following the diffuse domain equation derivation provided in Li et al. [18]. Similar derivation appeared in Ramière et al. [12]. To begin, multiply Eq. (1) by a test function ψ and integrate it over the fluid domain Ω f to obtain In the above equation we used the vector identity with the scalar field a = ψ and vector field b = κ ∇q, along with the Neumann boundary condition on the fluid-solid interface ∂Ω f as written in Eq.
This allows us to simplify the second integral in the left-hand side of Eq. (A.3) as Therefore, the weak form of the Poison equation in the extended domain can be written as Since ψ is an arbitrary test function, the collective term multiplying ψ in Eq. Next, noticing that ∇χ = δ ∂Ω f n, the forcing term of VP poisson equation becomes Substituting δ ∂Ω f g term from Eq. Defining a flux function β that satisfies the property of β · n = −g on ∂Ω f , the above equation can be written as The In order to avoid computing the gradient of a possible discontinuous indicator function χ, the above equation is re-written as (B.6) The normal vector appearing in the first term of Eq. (B.6) can be computed numerically using the signed distance function as n = −∇φ.

C. Order of accuracy results using Approach A with spatially varying flux values
In this section, we present the spatial convergence rate of the error norms using Approach A (analytical construction of β) for the egg and torus domains. The test problem remains the same as defined in Secs. 5.2 and 5.3 for the egg and torus domain, respectively. The flux-forcing function using Approach A is β = κ∇q exact . As shown in Fig C.13, we observe second-order convergence rate for both problems using Approach A.

D. Comparison of Approach C with Approach D
Here, we demonstrate the efficacy of Approach C in comparison to Approach D for numerically constructing the flux-forcing functions. To implement Approach D, we discretize the interface into a set of discrete Lagrangian/marker points with position X ≡ (X, Y ) and spread the two components (in 2D) of β(X) ≡ (β X , β Y ) to the nearby x-and y-faces of the Cartesian grid cells, respectively. We consider two kernel functions for the spreading operator: a one-point top hat function and a six-point Gaussian bell-like Here, x is the face-center location of the x-face and h is the grid cell size. Similarly, the one-dimensional form of the six-point spline function reads as (D. 2) in which σ = |r| + 3. The graphical representation of these functions is shown in Fig. D.14. In dimensions higher than one, a tensor-product form of the one-dimensional functions is used. We refer the readers to Peskin [1] for more details on the spreading operator. The distance between the maker points is kept approximately equal to h, although increasing or decreasing the distance up to a factor of two did not affect the overall accuracy of the scheme (data not shown). We consider a test problem similar to the one defined in Sec. 5.2, in which a circle of radius 3/2 is embedded into a larger computational domain of extents Ω ∈ [0, 2π] 2 . Inhomogeneous Neumann boundary conditions with g(x) = −κ n·∇q exact are imposed on the fluid-solid interface ∂Ω s , whereas Dirichlet boundary conditions are imposed on the external boundaries of the computational domain, i.e., q| ∂Ω(x) = q exact (x). We present the contours of the numerical (red) and exact (blue) solutions and the convergence rate of the numerical schemes based on Approach C and D in Fig. D.15. As observed in the figure, excellent agreement is obtained between the numerical and exact solutions in the fluid domain (considered to be outside the cylinder) with Approach C. However, when Approach D is used considering either the top hat or the spline function, there is a large disagreement between the numerical and exact solutions. Furthermore, Approach C exhibits approximately first-order accuracy using the continuous indicator function, whereas Approach D exhibits zeroth-order accuracy. A similar discrepancy is observed using the discontinuous indicator function with Approach D, although it is slightly less severe than the continuous case. In contrast, the order of accuracy of Approach C improves further when the discontinuous χ is used. Data for the discontinuous function is not shown in the interest of brevity. Based on the results of this section we do not recommend Approach D to impose the spatially varying Neumann/Robin boundary conditions.

E. Effect of the penalization parameter
In this section, we study the effect of the penalization parameter η on the order of accuracy of the fluxbased VP method. The test problem described in Sec. 5.2 for the hexagram interface is considered here. We solve this problem using four different η values: η = {10 −2 , 10 −4 , 10 −8 , 10 −12 }. As noted in Fig. E.16, except for the largest value of η = 10 −2 , the convergence rate remains the same for the rest of the η values. Based on the results of this section we chose η = 10 −8 for all our test cases.
F. Error norm and curve-fitting of the data