A diffuse interface method for simulation-based screening of heat transfer processes with complex geometries
AA diffuse interface method for simulation-basedscreening of heat transfer processes with complexgeometries
Elizabeth J. Monte a , James Lowman a , Nasser Mohieddin Abukhdeir a,b, ∗ a Department of Chemical Engineering b Department of Physics and AstronomyUniversity of Waterloo200 University Avenue WestWaterloo, Ontario, Canada N2L 3G1
Abstract
Frequently, the design of physicochemical processes requires screening of largenumbers of alternative designs with complex geometries. These geometries mayresult in conformal meshes which introduce stability issues, significant computa-tional complexity, and require user-interaction for their creation. In this work, amethod for simulation of heat transfer using the diffuse interface method to cap-ture complex geometry is presented as an alternative to a conformal meshing,with analysis and comparisons given. The methods presented include automatednon-iterative generation of phase fields from CAD geometries and an extensionof the diffuse interface method for mixed boundary conditions. Simple mea-sures of diffuse interface quality are presented and found provide predictions ofperformance. The method is applied to a realistic heat transfer problem andcompared to the traditional conformal mesh approach. It is found to enablereasonable accuracy at an order-of-magnitude reduction in simulation time orcomparable accuracy for equivalent simulation times.
Keywords: diffuse interface; design screening; complex domain; phase field;heat transfer analysis
1. Introduction
Computational multiphysics simulations are a key enabler for the advance-ment of the design and optimization of physicochemical processes. Design, con-trol, and optimization of these processes requires detailed insight, and low-cost nonhazardous methods for the evaluation of candidate designs, for whichsimulation-based screening is ideal. However, several challenges posed by mul-tiphysics simulations of industrially relevant processes, scales, and geometries ∗ Corresponding author, [email protected]
URL: https://uwaterloo.ca/comphys (Nasser Mohieddin Abukhdeir)
Preprint submitted to Elsevier January 19, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] J a n xist such as computational complexity and numerical stability (Slotnick et al.(2014)).Industrial processes involve geometries with both complex topology and sig-nificantly varying length scales. In order to perform computational multiphysicssimulations of these processes, standard spatial discretization involves the use ofunstructured conformal meshes , which capture domain shape through conform-ing to its boundary. This typically results in spatial discretization length scalesbeing far smaller than those characteristic to the underlying physicochemicalphenomena, but which are required to capture boundary curvature. Conse-quently, the numerical solutions of multiphysics models using conformal meshescan suffer from stability issues, eg. Courant number constraints for convection-dominated processes, and have significant additional computational complexitycompared to structured meshes of comparable size. Also, the computationalcost and frequent need for user interaction in the creation of conformal meshes(Franz et al. (2012)) are detrimental to the evaluation of large numbers of variedgeometries as part of a design process. This motivates the use of alternativemethods for mesh generation and, more generally, the imposition of complex do-main boundaries and boundary conditions (Anderson et al. (1998); Scardovelliand Zaleski (1999); Mittal and Iaccarino (2005)).There exist many approaches for approximating complex boundaries andinterfaces as an alternative to conformal meshing. In general, these methodsinvolve the creation of a nonconformal mesh and the imposition of boundaryconditions through the addition of various types of forcing terms within a subdo-main of the nonconformal mesh that conforms to the desired complex geometry.These methods have been classified in the literature depending on the type ofboundary or interface being modelled: solid/solid, fluid/solid or fluid/fluid.
Im-mersed boundary methods are those which were developed for solid/solid andfluid/solid interfaces (Mittal and Iaccarino (2005)).
Diffuse-interface methods are those which were developed for fluid/fluid interfaces (Anderson et al. (1998)).There is significant overlap and interchange of this classification and so-calledcontinuous-forcing immersed boundary methods (Mittal and Iaccarino (2005))and diffuse-interface methods are similar in approach.Over the past decade, substantial progress has been made with respect tothe accuracy, stability, and application of these types of methods (Nguyen et al.(2017)). Rami`ere et al. (2007) developed a continuous immersed boundarymethod for mixed boundary conditions (Dirichlet, Neumann, Robin). Li et al.(2009b) developed different continuous immersed boundary condition approx-imations (Dirichlet, Neumann, and Robin) compatible with a moving diffuseinterface, and showed that the resulting solutions converge to the conformalboundary condition in the sharp interface limit. This was accomplished throughthe use of the method of matched asymptotic expansions. Lerv˚ag and Lowen-grub (2015) performed an asymptotic analysis of Neumann and Robin boundaryconditions on a steady reaction-diffusion problem and compared the accuracyof different diffuse interface formulations. Results showed that the boundaryconditions proposed in Li et al. (2009b) are either first order or second orderaccurate with respect to diffuse interface length scale, depending on which for-2ulation is used. Schlottbom (2016) examined the diffuse interface on ellipticproblems with Dirichlet boundary conditions, showing rates of convergence tothe conformal solution for uniform and adaptive mesh refinement (local to theinterface). Franz et al. (2012) also showed the convergence of the Dirichletboundary conditions formulated by Li et al. (2009b), similar to the results pre-sented in Lerv˚ag and Lowengrub (2015), and again applied these methods tothe solution of a steady reaction-diffusion equation.The stability of diffuse interface methods has been advanced in the pastdecade, particularly the stability of Dirichlet boundary conditions (Juntunenand Stenberg (2009); Nguyen et al. (2018)), through the use of the Nitschemethod (Nitsche (1971)). The Nitsche method provides a stable variationallyconsistent method (Nguyen et al. (2018); Evans and Hughes (2013)) for weaklyenforcing Dirichlet boundary conditions by incorporating flux terms into theweak form of the differential equation to penalize deviations from those bound-ary conditions. Nguyen et al. (2018) show that a single parameter introducedin the Nitsche method should be chosen as the largest eigenvalue for that sys-tem. Schillinger et al. (2016) have recently developed a non-symmetric Nitschemethod, which is parameter-free, which avoids potential stability issues resultingfrom using extreme values for the stabilization parameter.The diffuse interface method has been applied to a broad range of appli-cations. Teigen et al. (2009) applied the method to modelling diffusion andadsorption/desorption within a deformable interface using both adaptive meshrefinement and a multigrid method. Aland et al. (2010) used the diffuse interfacemethod with the Cahn-Hilliard/Navier-Stokes model for immiscible two-phaseflow in complex geometries, accounting for contact lines. Nguyen et al. (2017)used the diffuse interface method to perform stress analysis of bone structuresdirectly from tomography data of bone mineral density. They formulated thediffuse interface through artificial diffusion of three-dimensional tomographydata and compared their results to standard voxel finite-cell methods, findingsimilar accuracy. Stoter et al. (2017) applied this approach to the simulation ofsingle-phase/porous media transport within the human liver using MRI imaging.However, they found that for equivalent accuracy, using a traditional conformalmesh approach was less computationally complex. Recently, Treeratanaphitakand Abukhdeir (2021) applied the diffuse interface method to the solution of thetwo-fluid model for liquid/dispersed gas flows, finding that the kernels used forconstructing the diffuse-interface had negligible effect on accuracy for interfacescales significantly smaller than the geometry scale.While significant progress has been made, there are still significant prac-tical challenges to the application of diffuse interface methods for engineeringdesign. Firstly, computation of a phase field corresponding to a geometry spec-ified through CAD software is a non-trivial task. Methods have been developedfor this task (Nguyen et al. (2017, 2018); Teigen et al. (2009); Aland et al.(2010); Stoter et al. (2017); Aland et al. (2020)), but all involve the solutionof a diffusion equation to generate a smooth differentiable phase field or othercomputationally complex algorithms. Additionally, physicochemical processesrarely involve a single uniform boundary condition and, instead, typically involve3arying boundary conditions over sub-domains of the conformal surface. Yet,most applications of diffuse interface methods have been limited to the uniformboundary condition case. Li et al. (2009b) presented mathematical formula-tions for Dirichlet, Neumann, and Robin boundary conditions, but only applieda single boundary condition in each of their test cases. Lerv˚ag and Lowengrub(2015) tested spatially varying Neumann and Robin boundary conditions, butdid not combine different types of boundary conditions in the same simulation.Rami`ere et al. (2007) combined Dirichlet, Neumann, and Robin boundary con-ditions in their test cases. However, the Neumann boundary conditions wereapplied to portions of the mesh which conform to the complex geometry bound-ary, not to areas approximated with a diffuse interface. They did, in later work,combine Dirichlet and Robin boundary conditions over different sections of adiffuse interface boundary approximation, although it is unclear how the specificboundary sections were isolated (Rami`ere et al. (2007)).Within the aforementioned context, the overall objective of this work isto develop a comprehensive diffuse interface method for use in screening heattransfer processes with complex geometries. That is, the intent is not to replaceconformal mesh simulations but instead to reduce the computational complex-ity during the screening process through approximation of the geometries andboundary conditions. Specific objectives include:1. Development of an automated and non-iterative method for the generationof phase fields from CAD descriptions of geometries, including segmentingdifferent boundary conditions.2. Development of a stable method for imposing mixed boundary conditionsusing the diffuse interface method.3. Validation of the developed methods with the solution of a heat transferproblem with spatially-varying boundary conditions on a complex geome-try and evaluate performance compared to the reference/conformal meshsolution.The paper is organized as follows: First, the diffuse interface and Nitschemethods are reviewed and demonstrated for Poisson’s equation in Section 2.Then, the method for generation of phase fields from CAD geometries andquality measures for predicting its performance are presented and demonstratedin Sections 3.1-3.2. A method for imposing mixed boundary conditions for thediffuse interface is then presented in Section 3.3 and applied to a heat transferproblem (LED heat sink) in Section 3.4. Finally, conclusions are summarizedin Section 4.
2. Model Formulation
A large and diverse amount of past research has been devoted to phasefield methods. Research on and applications of phase field methods has mainlyfocused either on fluid-structure interactions (Mittal and Iaccarino (2005)), orfluid-fluid interactions near and far from critical points (Anderson et al. (1998)).4ome previous work has also examined phase field methods for use on a singledomain, particularly Nguyen et al. (2018). In this work, the use of a stationarydiffuse interface is focused on and this section is limited to relevant background.
For a given problem, some complex geometry Ω is defined as the simulationdomain. The standard conformal meshing approach would be to discretize Ωwith an unstructured mesh which conforms to its boundary ( ∂ Ω), shown inFigure 1. As can be seen, this approach to capturing domain geometry does not,in general, exactly conform apart from boundary nodes of the mesh elements.Conformance increases as spatial discretization scale decreases and converges tothe exact geometry as this scale approaches zero. Additionally, generation ofconformal meshes typically requires them to be unstructured, involves significantcomputational complexity, and significantly increases the complexity of meshpartitioning for parallel computation. (a) (b) (c) (d)
Figure 1: An example complex geometry meshed by conformal meshes with increasing levelsof refinement. A high level of refinement is needed for mesh edges to reasonably conform tothe geometry boundary, particularly in regions of high curvature.
Alternatively, the diffuse interface method defines a domain κ , typicallyCartesian (rectangular cuboid), which encompasses Ω and is then discretizedwithout conforming to the boundaries ∂ Ω, as shown in Figure 2a. A structuredmesh may be used, which requires significantly less computational complexitycompared to unstructured meshes required by the conformal approach. In orderto capture the desired original geometry, a scalar phase field φ is defined on κ .The phase field has the value of one within Ω, the value zero outside of Ω,and may be discontinuous (sharp-interface) or continuous (diffuse interface) atthe interface between the two subdomains, shown in Figures 2b-c, respectively.The sharp boundary of the original geometry can be approximated by the levelset φ = 0 .
5, while the diffuse approximation extends over the region |∇ φ | > δ as the diffuseinterface width approaches zero, the sharp-interface limit, as shown in Figures2e-f. 5 a) (b) (c) (d)(e) (f) Figure 2: An example complex geometry approximated by the diffuse interface method with(a) the geometry encompassed by a structured mesh (mesh coarsened for visibility), (b) thephase field approximating the geometry on a refined structured mesh, (c) the same phase fieldnow with a diffuse interface, (d) the diffuse approximation of the interface indicated by |∇ φ | ,and further schematics of (e) the phase field and (f) its gradient for different diffuse interfacewidths. Based on Figures 1, 3, and 4 in Nguyen et al. (2018). Given a diffuse interface approximation of a domain Ω, the weak form ofthe well-posed problem may be reformulated from a conformal to diffuse inter-face (nonconformal) form using the following integral identities (Nguyen et al.(2018)), (cid:90) Ω A d
Ω = (cid:90) κ HA dκ ≈ (cid:90) κ φA dκ (1) (cid:90) ∂ Ω B d
Ω = (cid:90) κ δ ∂ Ω B dκ ≈ (cid:90) κ | ∇ φ | B dκ (2)where the phase field and its gradient approximate the Heaviside function H and the Dirac delta function δ ∂ Ω respectively. Functions A and B can representscalar or vector fields. The unit normal vector to the complex domain boundarycan also be approximated as follows (Nguyen et al. (2018)), n ≈ − ∇ φ | ∇ φ | (3)This approach is demonstrated using the Poisson problem with generalized6oundary conditions and formulated for the finite element method, −∇ u = f on Ω (4) u = h on ∂ Ω D (5) − n · ∇ u = g on ∂ Ω N (6) − n · ∇ u = r ( u − q ) on ∂ Ω R (7)where ∂ Ω D , ∂ Ω N , and ∂ Ω R denote the boundary sections on which Dirichlet,Neumann, and Robin boundary conditions are applied respectively. By intro-ducing a test (or weighting) function v in the Sobolev space H (Ω), the weakform can be formulated, (cid:90) Ω ( ∇ u + f ) v d Ω = 0 (8) (cid:90) Ω ∇ u · ∇ v d Ω − (cid:90) Ω ∇ · ( v ∇ u ) d Ω = (cid:90) Ω f v d Ω (9) (cid:90) Ω ∇ u · ∇ v d Ω − (cid:90) ∂ Ω v ( n · ∇ u ) d∂ Ω = (cid:90) Ω f v d Ω (10)where integration by parts and the divergence theorem are used.Neumann and Robin boundary conditions are imposed through substitutioninto the surface integral term, (cid:90) Ω ∇ u · ∇ v d Ω + (cid:90) ∂ Ω N vg d∂ Ω = (cid:90) Ω f v d Ω (11) (cid:90) Ω ∇ u · ∇ v d Ω + (cid:90) ∂ Ω R vr ( u − q ) d∂ Ω = (cid:90) Ω f v d Ω (12)Dirichlet boundary conditions are imposed by imposing both u = h and v = 0on ∂ Ω D .The diffuse interface method transforms the volume integrals on Ω and thesurface integrals on ∂ Ω into volume integrals on κ . Eqns. 1-3 can be insertedinto eqns. 11-12 to produce the diffuse interface formulation for Neumann andRobin boundary conditions, (cid:90) κ ∇ u · ∇ vφ dκ + (cid:90) κ vg | ∇ φ | dκ = (cid:90) κ f vφ dκ (13) (cid:90) κ ∇ u · ∇ vφ dκ + (cid:90) κ vr ( u − q ) | ∇ φ | dκ = (cid:90) κ f vφ dκ (14)Dirichlet boundary conditions require additional treatment, as will be discussedin the following section.Since φ goes to zero outside of the approximation to the complex geometry,singularities are likely if κ is much larger than Ω. This can be prevented byredefining the phase field as follows, φ (cid:48) = α + (1 − α ) φ (15)7here φ (cid:48) is the new phase field and α is some small constant eg. × − .Alternatively, any mesh elements in κ on which φ is zero can be removed beforesolving (Nguyen et al. (2018)). The diffuse interface method requires additional treatment of Dirichlet bound-ary conditions due to the lack of an explicit conformal boundary. Most previouswork on phase field methods has focused on penalty methods or variants of theNitsche method (Nitsche (1971)).The penalty method directly penalizes deviations from the Dirichlet bound-ary condition at the diffuse boundary by adding the term (Juntunen and Sten-berg (2009)), β (cid:90) ∂ Ω D v ( u − h ) d∂ Ω D (16)to the mathematical formulation, where β is a stabilization parameter whichis problem dependent. Li et al. (2009a) show several variants of the penaltymethod; Schlottbom (2016) and Rami`ere et al. (2007) also use penalty-typemethods. However, the Nitsche method (Nitsche (1971)) has been recentlyshown to enable higher accuracy and less sensitivity to the choice of a stabiliza-tion parameter compared to the penalty method (Nguyen et al. (2018)). TheNitsche method both penalizes deviation from the Dirichlet boundary conditionand includes additional terms to maintain symmetry (Benk et al. (2012)), (cid:90) Ω ∇ u · ∇ v d Ω − (cid:90) ∂ Ω D u ( n · ∇ v ) d∂ Ω D − (cid:90) ∂ Ω D v ( n · ∇ u ) d∂ Ω D + β (cid:90) ∂ Ω D v ( u − h ) d∂ Ω D = (cid:90) Ω f v d Ω − (cid:90) ∂ Ω D h ( n · ∇ v ) d∂ Ω D (17)where β is a stabilization parameter which is problem dependent.The Nitsche method was first applied within the phase field method byFreund and Stenberg (1995) and has since been applied in other studies (Nguyenet al. (2017, 2018); Stoter et al. (2017); Benk et al. (2012)). Using eqns. 1-3,the diffuse interface formulation incorporating the Nitsche method becomes, (cid:90) κ ∇ u · ∇ vφ dκ + (cid:90) κ u ( ∇ φ · ∇ v ) dκ + (cid:90) κ v ( ∇ φ · ∇ u ) dκ + β (cid:90) κ v ( u − h ) | ∇ φ | dκ = (cid:90) κ f vφ dκ + (cid:90) κ h ( ∇ φ · ∇ v ) dκ (18)The Nitsche method has also been extended to Neumann and Robin boundaryconditions by Juntunen and Stenberg (2009).As shown in eqns. 17-18, the Nitsche method introduces a new parameter β , generally defined as β = γ/h (Juntunen and Stenberg (2009); Evans andHughes (2013); Benk et al. (2012); Hansbo (2005); Ruess et al. (2013, 2014);de Prenter et al. (2018)). The characteristic length of each mesh element h , eg. the average or maximum scale of the mesh element. With a structured mesh,8uch as those used by the diffuse interface method, h is constant over the entiremesh, but an unstructured mesh would require h to be defined element-wise(de Prenter et al. (2018)). A stability constant γ , must be large enough tosatisfy a stability inequality, see Hansbo (2005) for further discussion. Evansand Hughes (2013) show that γ should be proportional to the square of thenumerical solution interpolant order. Benk et al. (2012), Ruess et al. (2013),and Ruess et al. (2014), further find that γ values of 2 −
10, 10 − − β , β = 10 n ∆ x (19)where n refers to the interpolant order of the numerical solution and the meshspacing ∆ x is used for h . β could also be calculated element-wise as the max-imum eigenvalue of the stability inequality, see, for example, (Nguyen et al.(2018); de Prenter et al. (2018)).
3. Results and Discussion
The two complementary methods developed in this work, within the contextof enabling the use of the diffuse interface method for screening of heat transferprocesses, are:1. An automated and non-iterative method for generating phase fields fromCAD geometries, including segmenting different boundary conditions.2. A method for imposing mixed boundary conditions using phase fields.These methods will be presented and validated using simple variations of steady-state diffusion equations with uniform and mixed Dirichlet and Neumann bound-ary conditions.The first validation case is defined on a circular domain ( ˜ R = 1) and waschosen due to the smoothness of the geometry boundary, − ∇ T = (cid:18) πα (cid:19) r sin (cid:18) πr α (cid:19) − πα cos (cid:18) πr α (cid:19) (20)with an exact solution, T = sin (cid:18) πr α (cid:19) (21)for compatible Dirichlet and Neumann boundary conditions. The parameter α is set to 0.05.The second validation case is defined on a rectangular domain ( ˜ L = ˜ W = 1), − ∇ T = 8 π sin (2 πx ) cos (2 πy ) (22)with an exact solution, T = sin (2 πx ) cos (2 πy ) (23)9or compatible Dirichlet and Neumann boundary conditions. This rectangulardomain has sharp corners, which introduces more challenging geometric condi-tions for the imposition of a (continuous) diffuse interface approximation to theboundary.Finally, the methods developed in this work were applied to a realistic heattransfer application, evaluation of the LED heat sink design shown in Fig-ure 3. This is a diffusion problem with mixed Neumann (see Figure 3b) andRobin boundary conditions, solved on a complex three dimensional geometry. Itdemonstrates the performance benefits of the diffuse interface method comparedto a conformal mesh solution. (a) (b) Figure 3: Example of a heat sink with a complex shape with (a) top-view in contact withthe environment and (b) bottom-view in contact with a spatially varying heat source (in thiscase, several LEDs). CAD is based on Ghosh (2015) with modifications.
All simulations detailed in the following sections used the Netgen/NGSolve(Sch¨oberl (2020)) implementation of the finite element method. Conformalmeshing of the complex geometries used Netgen/NGSolve algorithms for au-tomated mesh generation and defect removal. Algorithms to generate and markthe phase field were implemented in Python and have been provided as supple-mental material.
Past research applying phase field methods to scientific and engineering ap-plications (Nguyen et al. (2017); Rami`ere et al. (2007); Rami`ere et al. (2007);Nguyen et al. (2018); Teigen et al. (2009); Aland et al. (2010); Stoter et al.(2017); Aland et al. (2020); Benk et al. (2012); De (2016); Mo et al. (2018))has leveraged the ease of capturing complex, potentially deforming, geometrieswhere a conformal mesh would be computationally infeasible. In general, thesepast methods involve the projection of a binary discontinuous phase field, cre-ated from the geometry, onto the computational grid or mesh. The binaryphase value is then used to distinguish the inside/outside of the geometry as asharp-interface. However, this results in a discontinuity at the interface, shown10n Figure 2b, which is suitable for the application of discontinuous immersedboundary methods (Mittal and Iaccarino (2005)), but not for continuous im-mersed boundary methods.In general, a continuous phase field can be generated from the discontinu-ous phase field (see Figures 2c-d) through the solution of a diffusion equationor the use of a signed distance function.
Ad hoc “blurring” has been used inimmersed boundary methods to spread and interpolate the forcing functions(for example, Toja-Silva et al. (2014)) and could also be applied to the diffuseinterface method. However, the continuous phase fields produced by blurringmethods do not have easily quantifiable diffuse interface thicknesses, nor a uni-form thickness throughout the boundary, so it becomes challenging to enforcethe quality of the interface approximation. Diffusion equation approaches solvethe Allen-Cahn equation over the domain, using the discontinuous phase field asthe initial condition and specifying an interface length scale parameter to con-trol the width of the resulting phase field (for example, Nguyen et al. (2018)).This phase field varies continuously from zero to one, behaving as the hyperbolictangent function within the interface region. However, these methods requireiterating through a series of pseudo-timesteps to reach the metastable diffusedstate and additional computational complexity is required for numerical differen-tiation (assuming an explicit integration method is used). Since the metastablesolution of the Allen-Cahn equation is known as a function of the distance fromthe initial discontinuous interface, this iterative procedure could be avoided ifsuch a signed distance function were known. Several works have suggested theuse of a signed distance function (Franz et al. (2012); Li et al. (2009b); Lerv˚agand Lowengrub (2015); Schlottbom (2016); Teigen et al. (2009)), but have notsuggested a method for generating this signed distance function when it is notknown analytically. Aland et al. (2020) provide a method for computing thesigned distance function of any geometry from its surface triangulation. How-ever, their use of unstructured meshes necessitates the use of costly ray-tracingalgorithms.Similar to the work of Aland et al. (2020), this work provides a general non-iterative method suitable for any geometry. It leverages the diffuse interfacemethod’s structured meshes to use far more computationally efficient distancetransform algorithms. The algorithm is based on the combination of a distancetransform and a kernel function . The distance transform approximates a scalardistance function (from the interface) d ( x ). It is then used to generate thecontinuous phase field by applying the user-specified kernel function φ ( d ) thatdefines the structure of the continuous interface in closed form. The width ofthe continuous interface is controlled by a user-specified parameter λ . Pastwork (Treeratanaphitak and Abukhdeir (2021)), analyzed different options forthe kernel function and found minimal differences in solution accuracy when theinterface width was small with respect to domain size. The results presentedin this work use the basis/interpolant function itself as the kernel function, andthe error function was also considered. The full algorithm is summarized asfollows: 11. Given the geometry specified as a discretized surface (boundary faces,edges, and vertices), generate a structured quadrilateral mesh that fullyencompasses the geometry. This is the nonconformal mesh.2. Define a temporary phase field η ( x ) on the nonconformal mesh.3. Iterate through the geometry faces and set η to one at any mesh nodeswithin a reasonable tolerance of the face. This maps the border of thegeometry over the nonconformal domain.4. Calculate the chessboard distance transform (Fabbri et al. (2008)) of everymesh point in the nonconformal domain from this border, d ( x , η ).5. Scale the distance transform by the user-specified diffuse interface scale λ .6. Apply the user-specified kernel function to the scaled distance transformto produce a smooth transition from zero to one to zero over the borderregion. λ sets the width of this transition.7. Use a floodfill algorithm (Russ (2015)) to set η to one in the interior ofthe border generated in step 3. This maps a binary representation of thegeometry over the non-conformal domain.8. Combine the kernel function output and η to produce an phase field thatvaries from − φ ( x ).There are two parameters inherent to the diffuse interface approximation toa domain boundary. A third parameter is required when the Nitsche method isused to impose Dirichlet boundary conditions.Mesh Element Scale ∆ x : This can be calculated for a given mesh dimensionfrom the size of the nonconformal domain ( L ) and the number of meshelements ( N ) along said dimension, ∆ x = L/ ( N − x can be differentfor each different dimension of the nonconformal mesh. As shown in Figure4a, decreasing ∆ x increases the accuracy of the solution by increasing theorder of the spatial approximation.Diffuse Interface Scale λ : This is a measure of the diffuseness of the interface.The width of the diffuse interface, defined as the number of mesh elementsrequired for φ to change from zero to one, can be calculated as w = 2 λ/ ∆ x .As shown in Figure 4b, λ affects the accuracy of the solution. Largevalues of λ extend the boundary condition constraints well into the interiorof the complex geometry and generally decrease solution accuracy. Thehighest accuracy is typically observed when the diffuse interface is thesame width as the mesh element scale, as this is the closest approximationto a sharp interface. For the same reasons, solution accuracy also improvesas the width of the diffuse interface decreases relative to the scale of thenonconformal domain (see Figure 4c). However, when applying Neumannor Robin boundary conditions it may be preferable to use an interfacewidth of several mesh elements as this smooths the gradient of φ at theboundaries and reduces numerical approximation errors.12itsche Stabilization Parameter β : This parameter is required when using theNitsche method to impose Dirichlet boundary conditions. As discussed inSection 2.2, β can be obtained by solving an eigenvalue problem (Nguyenet al. (2018)) or calculated as β = 10 n / ∆ x , where n is the order ofinterpolant used (Evans and Hughes (2013); Benk et al. (2012); Ruesset al. (2013, 2014)). (a) (b) (c) Figure 4: The effects of (a) ∆ x , (b) the width of the diffuse interface, and (c) the width ofthe diffuse interface relative to the size of the domain on the accuracy of the diffuse interfacemethod. The non-iterative method was applied to validation eqn. 20. The exact solu-tion and its derivative (Figure 5a) were used to formulate uniform Dirichlet andNeumann boundary conditions respectively, and solutions were obtained from astandard conformal mesh and the diffuse interface method. The convergence be-haviour, global error versus number of mesh elements, of both methods is shownin Figures 5b-5c. In both cases, the conformal mesh reached a more accuratesolution than the diffuse interface method. In the case of uniform Neumannboundary conditions, the diffuse interface method also had a slightly lower con-vergence rate and converged at a smaller mesh size. This is expected, since thediffuse interface method uses an inherently poorer interface approximation thana conformal mesh.The diffuse interface method is expected to show significant performancebenefits in the computation time needed to achieve a given level of accuracy.However, for this simple two-dimensional circular geometry the computationaloverhead associated with the solver was sufficiently high to obscure any differ-ences between the conformal mesh solution and the diffuse interface solution.A more complex three dimensional problem is used in Section 3.4 in order toevaluate this aspect of the performance of the method.13 a) (b) (c)
Figure 5: (a) surface plot of the exact solution (eqn. 21). Convergence behaviour of theconformal solution and diffuse interface method with (b) Dirichlet boundary conditions and(c) Neumann boundary conditions.
In the previous validation case, the geometry involved was a simple circlewith a constant and bounded mean curvature, defined as, κ = − ∇ s · n (24)where n is the unit normal to the boundary surface and ∇ s is the surfacedivergence. The second validation case involves a square geometry, which ischosen due to the presence of points on the bounding surface in which the meancurvature is singular, that is, there are corners in the geometry. This presents amore challenging geometry for the application of the diffuse interface method.The method was applied to validation eqn. 22, again for uniform Dirichletand Neumann boundary conditions obtained from the exact solution and itsderivative respectively (Figure 6a). The convergence behaviour, global errorversus number of mesh elements, for the diffuse interface method and a stan-dard conformal mesh are shown in Figures 6b-6c. Again, the conformal meshsolution shows better convergence behaviour in both cases. Particularly for theuniform Neumann boundary conditions, it is significantly more accurate andhas a significantly higher convergence rate than the diffuse interface method.Moderate differences are expected, but clearly the presence of discontinuities inthe curvature of the boundary has a further significant effect on the performanceof the diffuse interface method. The performance of the conformal mesh solu-tion is also especially biased in this case, since the domain geometry (square) isconformal to the mesh geometry (rectangular).14 a) (b) (c) Figure 6: (a) surface plot of the exact solution (eqn. 23). Convergence behaviour of theconformal solution and diffuse interface method with (b) Dirichlet boundary conditions and(c) Neumann boundary conditions.
The effects of geometry curvature on the diffuse interface method are furtherevident from the spatial distribution of error. Figure 7 shows surface plots ofthe relative error for each of the two test cases. Dark areas indicate regionsof high error and the intensity is scaled by the maximum error in each figuredomain. In the case of the circle geometry, error is distributed relatively uni-formly around the entire circular subdomain. It does increase slightly in areaswhere the geometry curvature is more poorly approximated and in areas wherethe exact solution has significant spatial variation (outer region). The squaregeometry, likewise, shows higher error along the two boundaries that intersectthe greatest magnitude of the cosine wave. However, the error is distributednonuniformly with the majority of it located at the corners of the square, dueto the singularities in the boundary curvature.
The previously presented results show that there is significant variation inthe performance of the diffuse interface method depending on the attributes ofthe complex geometry. Consequently, a set of simple predictive scalar measuresare developed to evaluate the accuracy of the method and guide the selection ofuser-defined scales (mesh and diffuse interface). Based a simple scaling analysis,several dimensionless parameters are proposed, formulated to range from 0 →∞ , where values close to zero predict that the diffuse interface method willperform well, while values approaching ∞ predict the opposite. They serve aspredictive measures of the potential accuracy of the diffuse interface methodbased on domain geometry and physical scales:Interface Scale Measure λ d /λ i : This is the ratio of the diffuse interface scale,a user-specified parameter, to the physical scale of the interface. Apartfrom simulations under conditions approaching critical points Andersonet al. (1998), this value will approach ∞ for most continuum mechanicalmodels.Mean Chord Scale Measure λ d /λ g : This is the ratio of the diffuse interfacescale to the mean scale of the geometry.15 a) (b) Figure 7: Surface plots of the local error for the diffuse interface method applied to the (a)circlular geometry and (b) square geometry test cases. In both cases the error is normalizedto lie on [0,1].
Minimum Chord Scale Measure λ d /λ n : This is the ratio of the diffuse interfacescale to the smallest of the necks present in the geometry.Curvature Scale Measure λ d /λ c : This is the ratio of the diffuse interface scaleto the mean radius of curvature of the domain, excluding discontinuitiesin curvature (edges, points).Discontinuity Scale Measure (cid:80) i L i /λ d : This is the ratio of the length of thediscontinuities in the curvature of the domain (for an edge the length ofthe edge, for a point L = λ d ) to the diffuse interface scale.The various measures were calculated for the square and circle geometry testcases. Theoretical expressions for mean chord length and radius of curvatureare known for these simple geometries, but would not be known for a generalcomplex geometry. Thus, the chord lengths and radii of curvature are insteadapproximated using the boundary vertices of the CAD representations of eachgeometry. As can be seen in Table 1, the approximated measures agree wellwith the theoretical measures. The minor discrepancies are due to the CADnot perfectly approximating the geometry boundaries. The minimum chordscale measure is not given in either case, since neither domain contains necks.Similarly, the circle geometry contains no discontinuities, so the discontinuityscale measure is zero. In contrast, the square domain contains a discontinuityat each of the four corners. The remaining portions of the square boundary16re straight lines, resulting in a curvature scale measure of zero. The measuresare sensitive to the smoothness of the CAD representation. For example, asufficiently coarse CAD representation of a circle geometry could give a nonzerovalue for the discontinuity scale measure. This is a desirable property, since thediffuse interface approximation of the geometry boundary is entirely dependenton the fidelity of the CAD representation. Table 1: Quality metrics for the two sample geometries.
Circle SquareMeasure Theory Calculated Theory CalculatedMean Chord Scale 0 . λ . λ . λ . λ Minimum Chord Scale N/A N/A N/A N/ACurvature Scale 0 . λ . λ The majority of experimentally and industrially relevant design problems in-volve spatially varying boundary conditions. For example, standard models forpressure-driven flow through a conduit involve the imposition of boundary condi-tions on velocity which vary from no-slip conditions on the conduit walls (Dirich-let velocity and Neumann pressure condition) to nonzero pressure/velocity inletand outlet conditions. To date, the vast majority of applications of the diffuseinterface method are restricted to uniform boundary conditions, or account forspatial variation in a non-generalized manner (Nguyen et al. (2017); Rami`ereet al. (2007); Nguyen et al. (2018); Teigen et al. (2009); Aland et al. (2020); Benket al. (2012); De (2016); Mo et al. (2018)). Consequently, a method is developedfor the imposition of spatially-varying boundary conditions that is complemen-tary to the method presented in the previous section for the generation of thediffuse interface itself.Similar to the requirement that the bounding surface of the complex geome-try must be specified, so too must the sub-surfaces for each individual boundarycondition. Given a set of i sub-surfaces, the proposed method introduces a setof scalar fields { φ i } in addition to the boundary phase field φ . These are usedto indicate the local boundary condition in the following general form, (cid:88) i φφ i ( A i n i · ∇ i T + B i T + C i ) = 0 (25)The values of { A i , B i } specify the type of boundary condition. A i , B i (cid:54) = 0 isa Robin boundary condition, A i = 0 is a Dirichlet boundary condition, and B i = 0 is a Neumann boundary condition. Eqn. 25 replaces the set of bound-ary conditions from the original formulation and is incorporated into the weakformulation as described in Section 2.1. In order to maintain stability for theimposition of Dirichlet boundary conditions the Nitsche reformulation is also17sed (Nitsche (1971)), which reformulates the Dirichlet condition into a Robincondition as shown in Section 2.2.The non-iterative method proposed is decomposed into (i) generation of theset of scalar fields { φ i } and (ii) modification of the diffuse interface formu-lation to impose boundary conditions using this set. The algorithm for gen-eration of the scalar fields is summarized as follows and demonstrated for atwo-dimensional domain in Figure 8:1. Given the geometry specified as a discretized surface (boundary faces,edges, and vertices), identify sub-sets of the discretized surface associatedwith each of the sub-surfaces on which an individual boundary conditionis to be imposed.2. Calculate the centroid of the discretized geometry or use a user-specifiedlocation.3. Define a set of scalar fields { φ i } on the nonconformal mesh. In two dimen-sions, each field is one at every mesh node inside the solid angle extendingfrom the centroid through both bounding vertices to the edges of the non-conformal mesh. In three dimensions this is the cone extending from thecentroid through all boundary section boundary vertices to the edges ofthe nonconformal mesh. The fields are zero everywhere else, resulting inmasks over only specific regions of the domain.4. Optionally (if specified by the user), apply a distance transform weightedby λ overlap to the edge of each mask, then apply the kernel function to theresult. This results in a diffuse transition between neighbouring boundaryconditions. (a) (b) (c) (d) Figure 8: (a) φ , (b) |∇ φ | for the same complex geometry, (c) the mask used to apply aboundary condition to a section of |∇ φ | , and (d) the same mask if the different boundaryconditions are made to diffuse into each other. The method was applied to the same two-dimensional test cases used insection 3.1, with eqns. 20-22 formulated for mixed Dirichlet and Neumannboundary conditions. The placement of the different boundary conditions isshown in Figure 9a,b. Figure 9e,f then shows the convergence behaviour of thediffuse interface method with spatially-varying boundary conditions, comparedto solutions on standard conformal meshes. As expected, the conformal solu-tions show better convergence behaviour, higher accuracy, except on very coarse18eshes, and higher rates of convergence. Indeed, in both cases, the diffuse inter-face method shows similar trends and accuracy for spatially-varying boundaryconditions as for uniform conditions. This shows, qualitatively, that the addi-tions to the diffuse interface method for imposing spatially-varying boundaryconditions have no significant effect on the performance of the method.19 a) (b)(c) (d)(e) (f)
Figure 9: The placement of the different boundary conditions on (a) the circle geometryand (b) the square geometry. (c,d) exact solutions (eqns. 21 and 23) and (e,f) convergencebehaviour for the diffuse interface method with mixed Dirichlet and Neumann boundary con-ditions. .4. Application to a heat transfer problem The methods presented in Sections 3.1-3.3 were then used for a realistic ap-plication; evaluation of the heat transfer performance of a three dimensionalLED heat sink design. This demonstrates the potential performance of thediffuse interface method with mixed boundary conditions on a complex geom-etry and allows for a full comparison with standard conformal mesh solutions,including a timing comparison.A CAD representation of the design is shown in Figure 3. The heat sinkitself has a diameter of 12 cm and a height of 5 .
59 cm. It is assumed to be purealuminum, with a heat capacity of 910 Jkg − K − and a thermal conductivityof 205 Wm − K − (Young et al. (2008)). The bottom of the heat sink is incontact with three rectangular LEDs, each of which output 100 W and whichhave areas of 6 .
45 cm , 1 .
61 cm , and 9 .
68 cm . The remaining finned surfacesexperience natural convection in air, with a convective heat transfer coefficientof 10 Wm − K − and an environmental temperature of 298 K. Steady-state heattransfer within the heat sink is modelled using the following approximations:1. Temperature-independent material properties.2. Convective heat transfer boundary conditions for the heat sink/air inter-face (Figure 3a), n · q = h ( T − T ∞ ) (26)which corresponds to a Robin boundary condition.3. Heat flux boundary conditions for the LED/heat sink interface (Figure 3b), which corresponds to a Neumann boundary condition.The complexity of the CAD geometry can be assessed with the mesh qualitymetrics described in Section 3.2. The interface scale measure is ∞ since thephysical interface is assumed to be infinitesimally thin. The mean chord scalemeasure is λ/ .
358 cm and the minimum chord scale measure is λ/ .
941 cm.The curvature scale measure is 0 as the heat sink has many flat edges. Finally,the discontinuity scale measure is 253 . /λ and is only comprised of edgediscontinuities. Comparing these values to Table 1 confirms that the LED heatsink is far more difficult for the diffuse interface method to approximate thana circle or square. All of the diffuse interface solutions use a diffuse interfacethickness of 0 .
165 cm, so the mesh quality metrics are then ∞ , 0 .
461 cm − ,0 . − , 0, and 1537 . × − . As shown in Figure 10a, the standard conformal mesh approachgives higher accuracy than the diffuse interface method for a given number ofmesh elements. This is expected, as the diffuse interface method approximationof a geometry boundary is, in general, less accurate than that of a conformalmesh. However, as shown in Figure 10b, the total simulation time required bythe diffuse interface method scales significantly better than that of the stan-dard approach. This is predominantly due to the increased computation timerequired to generate an unstructured conformal mesh compared to a structuredgrid. The use of solvers optimized for structured meshes would further increasethis timing difference, but was not feasible for this work. Accounting for totalsimulation time (see Figure 10c), the diffuse interface method provides similaraccuracy for a lower computation time, until very fine meshes are required. Asshown in Figure 10d, the diffuse interface method also gives reasonable esti-mates of key design parameters (in this case the maximum temperature of theheat sink) significantly faster than the standard conformal solutions. This isdesirable for engineering design activities where screening of a large number ofcandidate designs is only feasible for a constrained amount of simulation time.22 a) (b)(c) (d) Figure 10: (a) convergence behaviour of the conformal solution and the diffuse interfacemethod, (b) timing comparison of the conformal solution and the diffuse interface method,(c) timing compared to error, and (d) timing compared to the maximum temperature of theheat sink.
4. Conclusions
A comprehensive method for simulation of heat transfer processes using thediffuse interface to capture complex domain boundaries was presented and com-pared to a traditional conformal mesh approach. This method was developedas an alternative to traditional conformal mesh simulations for screening of alarge number of design variations as part of a design process. It reduces thecomputational complexity and avoids the user-interaction required by confor-mal mesh creation. However, the reduced computational complexity enabled bythe diffuse interface method comes at the cost of accuracy.The methods presented and demonstrated in this work include automatednon-iterative generation of phase fields from CAD geometries (Section 3.1) and23n extension of the diffuse interface method to impose mixed boundary con-ditions (Section 3.3). Additionally, simple measures of diffuse interface qualitybased on CAD and mesh properties were presented and demonstrated to providean understanding of the performance of the diffuse interface method.The diffuse interface method was then applied to a realistic heat transferproblem, evaluation of the performance of an LED heat sink, and comparedto the traditional conformal mesh approach.. It is found to enable reasonableaccuracy at an order-of-magnitude reduction in simulation time or comparableaccuracy for equivalent simulation times. This supports its use for screeninglarge numbers of design variations as part of the design heat transfer processesand its future extension to fluid dynamics, mass transfer, and multiphysicsprocesses.
Acknowledgements
This research was supported by the Natural Sciences and Engineering Re-search Council (NSERC) of Canada, the Ontario Ministry of Ministry of Re-search, Innovation and Science, and Compute Canada.
References
Aland S, Lowengrub J, Voigt A. Two-phase flow in complex geometries: Adiffuse domain approach. Comput Model Eng Sci 2010;57:77–106. doi: .Aland S, Stenger F, Muller R, Deutsch A, Voigt A. A phase field approach totrabecular bone remodeling. Front Appl Math Stat 2020;6:12. URL: . doi: .Anderson DM, McFadden GB, Wheeler AA. Diffuse-interface methods in fluidmechanics. Annu Rev Fluid Mech 1998;30:139–65. doi: .Benk J, Ulbrich M, Mehl M. The Nitsche method of the Navier-Stokes equa-tions for immersed and moving boundaries. In: Proceedings of the SeventhInternational Conference on Computational Fluid Dynamics, ICCFD7. Inter-national Conference on Computational Fluid Dynamics. Big Island, Hawaii;2012. .De A. A diffuse interface immersed boundary method for convective heat andfluid flow. Int J Heat Mass Transfer 2016;92:957–69.Evans JA, Hughes TJR. Isogeometric divergence-conformaing B-splines forthe steady Navier-Stokes equations. Math Models Methods Appl Sci2013;23(8):1421–78. doi: .24abbri R, Costa LDF, Torelli JC, Bruno OM. 2d euclidean distance trans-form algorithms: A comparative study. ACM Computing Surveys 2008;40(1).doi: .Franz S, Roos HG, G¨artner R, Voigt A. A note on the convergence analysis ofa diffuse-domain approach. Comput Methods Appl Math 2012;12(2):153–67.Freund J, Stenberg R. On weakly imposed boundary conditions for second orderproblems. In: Proceedings of the Ninth Int. Conf. Finite Elements in Fluids.Venice; 1995. p. 327–36.Ghosh A. Brushed aluminium heat sink for circular cobled. Online; 2015. URL: https://grabcad.com/library/brushed-aluminium-heat-sink-for-circular-cob-led-1 .Hansbo P. Nitsche’s method for interface problems in computational mechanics.GAMM-Mitteilungen 2005;28(2):183–206. doi: https://doi.org/10.1002/gamm.201490018 .Juntunen M, Stenberg R. Nitsche’s method for general boundary conditions.Math Comp 2009;78:1353–74.Lerv˚ag KY, Lowengrub J. Analysis of the diffuse-domain method for solv-ing PDEs in complex geometries. Commun Math Sci 2015;13(6):1473–500. URL: https://dx.doi.org/10.4310/CMS.2015.v13.n6.a6 . doi: .Li X, Lowengrub J, R¨atz A, Voigt A. Solving PDEs in complex geometries: Adiffuse domain approach. Commun Math Sci 2009a;7:81–107. doi: .Li XG, Liu DX, Xu SM, Li H. CFD simulation of hydrodynamics of valve tray.Chem Eng Process Process Intensif 2009b;48:145–51. doi: .Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech2005;37:239–61. doi: .Mo H, Lien FS, Zhang F, Cronin DS. An immersed boundary method for solv-ing compressible flow with arbitrarily irregular and moving geometry. Int JNumer Methods Fluids 2018;88(5):239–63. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.4665 . doi: . arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.4665 .Nguyen L, Stoter S, Baum T, Kirschke J, Ruess M, Yosibash Z, SchillingerD. Phase-field boundary conditions for the voxel finite cell method:Surface-free stress analysis of ct-based bone structures. Int J Numer Meth-ods Biomed Eng 2017;33(12):e2880. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.2880 . doi: . arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.2880 ;e2880 cnm.2880. 25guyen LH, Stoter SK, Ruess M, Sanchez Uribe MA, Schillinger D. Thediffuse Nitsche method: Dirichlet constraints on phase-field boundaries. IntJ Numer Methods Eng 2018;113(4):601–33. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5628 . doi: . arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5628 .Nitsche J. ¨Uber ein variationsprinzip zur l¨osung von Dirichlet-problemenbei verwendung von teilr¨aumen, die keinen randbedingungen unterworfensind. Abhandlungen aus dem Mathematischen Seminar der Universit¨atHamburg 1971;36(1):9–15. URL: https://doi.org/10.1007/BF02995904 .doi: .de Prenter F, Lehrenfeld C, Massing A. A note on the stability parameter inNitsche’s method for unfitted boundary value problems. Comput Math withAppl 2018;75(12):4322–36. doi: .Rami`ere I, Angot P, Belliard M. A fictitious domain approach with spreadinterface for elliptic problems with general boundary conditions. ComputMethods Appl Mech Engrg 2007;196:766–81. doi: .Rami`ere I, Angot P, Belliard M. A general fictitious domain method withimmersed jumps and multilevel nested structured meshes. J Comput Phys2007;225(2):1347 –87. URL: . doi: https://doi.org/10.1016/j.jcp.2007.01.026 .Ruess M, Schillinger D, Bazilevs Y, Varduhn V, Rank E. Weakly enforcedessential boundary conditions for NURBS-embedded and trimmed NURBSgeometries on the basis of the finite cell method. Int J Numer Methods Eng2013;95(10):811–46. doi: .Ruess M, Schillinger D, Ozcan AI, Rank E. Weak coupling for isogeometric anal-ysis of non-matching and trimmed multi-patch geometries. Comput MethodsAppl Mech Eng 2014;269:46–71. doi: .Russ JC. The image processing handbook. volume 6. CRC press, 2015.Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and inter-facial flow. Annu Rev Fluid Mech 1999;31:567–603. doi: .Schillinger D, Harari I, Hsu MC, Kamensky D, Stoter SK, Yu Y,Zhao Y. The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in im-mersed finite elements. Comput Methods Appl Mech Eng 2016;309:625–52. URL: . doi: https://doi.org/10.1016/j.cma.2016.06.026 .26chlottbom M. Error analysis of a diffuse interface method for elliptic prob-lems with Dirichlet boundary conditions. Appl Numer Math 2016;109:109–22.doi: .Sch¨oberl J. Netgen/NGSolve. https://ngsolve.org/ ; 2020. URL: https://ngsolve.org/ ; v. 6.2.2009.Slotnick J, Khodadoust A, Alonso J, Darmofal D, Gropp W, Lurie E, MavriplisD. CFD vision 2030 study: a path to revolutionary computational aero-sciences. Technical Report; National Aeronautics and Space Administration(USA); 2014.Stoter SK, M¨uller P, Cicalese L, Tuveri M, Schillinger D, HughesTJ. A diffuse interface method for the Navier–Stokes/Darcy equa-tions: Perfusion profile for a patient-specific human liver basedon MRI scans. Comput Methods Appl Mech Eng 2017;321:70– 102. URL: . doi: https://doi.org/10.1016/j.cma.2017.04.002 .Teigen KE, Li X, Lowengrub J, Wang F, Voigt A. A diffuse-interface approachfor modeling transport, diffusion and adsorption/desorption of material quan-tities on a deformable interface. Commun Math Sci 2009;4(7):1009.Toja-Silva F, Favier J, Pinelli A. Radial basis function (RBF)-based interpo-lation and spreading for the immersed boundary method. Comput Fluids2014;105:66–75. doi: .Treeratanaphitak T, Abukhdeir NM. Phase-bounded finite elementmethod for two-fluid incompressible flow systems; 2021. doi:10.1016/j.ijmultiphaseflow.2019.04.024