A diffuse interface box method for elliptic problems
AA diffuse interface box method for elliptic problems
G. Negrini a , N. Parolini a and M. Verani a January 28, 2021 a MOX, Dipartimento di Matematica, Politecnico di Milano, Piazza Leonardo da Vinci 32, I-20133 Milano,Italy
Abstract
We introduce a diffuse interface box method (DIBM) for the numerical approximation on com-plex geometries of elliptic problems with Dirichlet boundary conditions. We derive a priori H and L error estimates highlighting the rˆole of the mesh discretization parameter and of the diffuseinterface width. Finally, we present a numerical result assessing the theoretical findings. Keywords : box method, diffuse interface, complex geometries
The finite volume method (FVM) is a popular numerical strategy for solving partial differential equa-tions modelling real life problems. One crucial and attractive property of FVM is that, by construction,many physical conservation laws possessed in a given application are naturally preserved. Besides,similar to the finite element method, the FVM can be used to deal with domains with complex ge-ometries. In this respect, one crucial issue is the construction of the computational grid. To facethis problem, one can basically resort to two different types of approaches. In the first approach, amesh is constructed on a sufficiently accurate approximation of the exact physical domain (see, e.g.,isoparametric finite elements [9], isogeometric analysis [10], or Arbitrary Lagrangian-Eulerian formula-tion [11,16,17]), while in the second approach (see, e.g., Immersed Boundary methods [19], the PenaltyMethods [2], the Fictitious Domain/Embedding Domain Methods [4–6], the cut element method [7, 8]and the Diffuse Interface Method [18]) one embeds the physical domain into a simpler computationalmesh whose elements can intersect the boundary of the given domain. Clearly, the mesh generationprocess is extremely simplified in the second approach, while the imposition of boundary conditionsrequires extra work. Among the methods sharing the second approach, in this paper we focus on thediffuse interface approach developed in [20]. In parallel, we consider, for its simplicity, the piecewiselinear FVM, or box method , that has been the object of an intense study in the literature (see, e.g.,the pioneering works [3, 15] and the more recent [12, 21]).The goal of this paper is to propose and analyse a diffuse interface variant of the box method, inthe sequel named DIBM (diffuse interface box method), obtaining a priori H and L error estimatesdepending both on the discretization parameter h (dictating the accuracy of the approximation of thePDE) and the width (cid:15) of the diffuse interface (dictating the accuracy of the domain approximation).Up to our knowledge, this is new in the literature. Besides, the study of DIBM for elliptic problems,despite its simplicity, opens the door to the study of more complicated differential problems and tothe analysis of diffuse interface variants of more sophisticated finite volume schemes.The outline of the paper is as follows. In section 2 we briefly recall the box method, while insection 3, we present the diffuse interface box method (DIBM) along with a priori error estimates.Finally in section 4 we will provide a numerical test to validate the theoretical results. The numericalresults have been obtained using the open-source library OpenFOAM ® .1 a r X i v : . [ m a t h . NA ] J a n The box method
In this section, we recall (see [3, 15, 21]) the box method for the solution of an elliptic problem. Let D ⊂ R be a polygonal bounded domain (in the following section this hypothesis will be relaxed). Weconsider the following problem: (cid:40) − ∆ u = f, in Du = g, on Γ = ∂D, (2.1)where f ∈ L (Ω) and g ∈ H / (Γ).Let T h = { t i } be a conforming and shape regular triangulation of D . We denote by h t the diameterof t ∈ T h and we introduce the set V h = { v i } of vertices of T h with V h = V ∂h ∪ V oh , the set V oh containingthe interior vertices of T h . We denote by w v the set of triangles sharing the vertex v . On T h we definethe space of linear finite elements V h,g h = (cid:8) v h ∈ C ( ¯ D ) : v h | t ∈ P ( t ) ∀ t ∈ T h and v h = g h on ∂D (cid:9) , where g h is a suitable piecewise linear approximation of g on ∂D .Let B h = { b v } v ∈ V oh be the “box mesh” (or dual mesh) associated to T h . Each box b v is a polygonwith a boundary consisting of two straight lines in each related triangle t ∈ w v . These lines are definedby the mid-points of the edges and the barycentres of the triangles in w v .On B h we introduce the space of piecewise constant functions, W h = (cid:8) w h ∈ L ( D ) : w h ∈ P ( b v ) ∀ b v ∈ B h (cid:9) . The box method for the approximation of (2.1) reads as follows: find u B,h ∈ V h,g h such that a T h ( u B,h , w h ) = ( f, w h ) D ∀ w h ∈ W h , (2.2)where a T h ( u B,h , w h ) = − (cid:88) v ∈ V oh (cid:90) ∂b v ∂ v h ∂ n b w h d s, (2.3)being n b the outer normal to b v and ( · , · ) D is the usual L scalar product on D .Note that there holds (see [3, 15] for the two dimensional case and [21] for the extension to anydimension) (cid:90) ∂b v ∂ φ v (cid:48) ∂ n b w h d s = (cid:90) D ∇ φ v · ∇ φ v (cid:48) d x, ∀ v ∈ V oh , ∀ v (cid:48) ∈ V h , (2.4)where φ v is the usual hat basis function with support equal to w v .The relation (2.4) is crucial to show the following perturbation results (see [15, 21]): (cid:107)∇ ( u B,h − u G,h ) (cid:107) L ( D ) ≤ Ch (cid:107) f (cid:107) L ( D ) , (cid:107) u B,h − u G,h (cid:107) L ( D ) ≤ Ch (cid:107) f (cid:107) L ( D ) , (2.5)where h = max t ∈T h h t and u G,h ∈ V h,g h is the linear finite element approximation to the solution ofproblem (2.1). The aim of this section is to introduce a variant of the box method for the approximate solution ofproblem (2.1) in case of a general (non-polygonal) domain D ⊂ R , where in the spirit of [20] theDirichlet boundary condition is treated with a diffuse interface approach. To this aim we introduce anhold-all domain Ω such that D ⊂ Ω. In the sequel we will work under the hypothesis Γ = ∂D ∈ C , .With a slight abuse of notation we denote by T h a shape regular triangulation of Ω. It is worth notingthat T h is not conforming with D . Following [20] we first select a tubular neighbourhood S (cid:15) of Γ,where (cid:15) denotes the width of S (cid:15) (see Figure 1). Then we introduce the set S (cid:15)h which contain s all the2igure 1: Diffuse interface representation: D is a surrogate domain of Ω; Γ is the Dirichlet boundaryand S (cid:15) is its tubular neighbour.triangles of T h having non-empty intersection with S (cid:15) . Note that the width of the discrete tubularneighbourhood S (cid:15)h is δ + (cid:15) where δ is the maximum diameter of triangles crossed by ∂S (cid:15) .To proceed, we assume that there exists an extension ˜ g ∈ H (Ω) of the boundary data g.We set D (cid:15)h = D \ S (cid:15)h and introduce the function u (cid:15),h ∈ H ( D (cid:15)h ) such that u (cid:15),h = g on ∂D (cid:15)h , whichsolves the following continuos problem: (cid:90) D (cid:15)h ∇ u (cid:15),h · ∇ v = (cid:90) D (cid:15)h f v ∀ v ∈ H ( D (cid:15)h ) . (3.1)The solution u (cid:15),h is then extended to S (cid:15)h by setting u (cid:15),h = ˜ g in S (cid:15)h .The following results have been proved in [20, Thm 1.2]:1 (cid:15) + δ (cid:13)(cid:13)(cid:13) u − u (cid:15),h (cid:13)(cid:13)(cid:13) L ( D ) + 1 √ (cid:15) + δ (cid:13)(cid:13)(cid:13) ∇ u − ∇ u (cid:15),h (cid:13)(cid:13)(cid:13) L ( D ) ≤ C (cid:16) (cid:107) f (cid:107) L ( D ) + (cid:107) g (cid:107) H ( D ) (cid:17) . (3.2)Let V (cid:15)h, ˜ g h = (cid:8) v h | D (cid:15)h : v h ∈ P ( t ) ∀ t ∈ T h and v h = ˜ g h on ∂D (cid:15)h (cid:9) , with ˜ g h the Lagrangian piecewiselinear interpolant of ˜ g .It has been proved (cf. [20, Thms 5.1 and 5.3]) that the linear finite element approximation u (cid:15)G,h ∈V (cid:15)h, ˜ g h of u (cid:15),h satisfies the following estimates: (cid:13)(cid:13)(cid:13) ∇ ( u (cid:15),h − u (cid:15)G,h ) (cid:13)(cid:13)(cid:13) L ( D ) ≤ C ( √ δ + κ + h ) (cid:16) (cid:107) f (cid:107) L ( D ) + (cid:107) ˜ g (cid:107) H ( D ) (cid:17) , (cid:13)(cid:13)(cid:13) u (cid:15),h − u (cid:15)G,h (cid:13)(cid:13)(cid:13) L ( D ) ≤ C ( δ + κ + h ) (cid:16) (cid:107) f (cid:107) L ( D ) + (cid:107) ˜ g (cid:107) H ( D ) (cid:17) , (3.3)where κ is the maximum diameter of the triangles intersection ∂S (cid:15) + h and u (cid:15)G,h has been extended to D (cid:15)h by setting u (cid:15)G,h = ˜ g h on S (cid:15)h .Let us now introduce the box method with diffuse interface (DIBM). We denote by u (cid:15)B,h ∈ V (cid:15)h, ˜ g h ,the approximation obtained from applying the box method to (3.1) (cf. (2.2)). The solution u (cid:15)B,h isthen extended to D by setting u (cid:15)B,h = ˜ g h in S (cid:15)h . Then employing the triangle inequality in combinationwith (3.2), (3.3) and (2.5) we get the following estimates for DIBM: (cid:13)(cid:13) ∇ ( u − u (cid:15)B,h ) (cid:13)(cid:13) L ( D ) (cid:46) √ (cid:15) + δ + √ δ + k + h, (cid:13)(cid:13) u − u (cid:15)B,h (cid:13)(cid:13) L ( D ) (cid:46) (cid:15) + δ + k + h . (3.4) In this section we numerically assess the theoretical estimates obtained in Section 3. To this aim,we consider the test case originally introduced in [20, Section 6] that is briefly recalled in the sequel.3igure 2: Discrete diffuse interface representation on triangulation (left) and on box mesh (right).Constrained cells are marked with red dots while the continuous and discrete diffuse interfaces arecoloured by darker an lighter red respectively.Let Ω = ( − , and let Γ be the boundary of the circle B (0) with centre (0 ,
0) and unitary radius.Thus, Γ splits the domain Ω into two subregions: D = B (0) and D = Ω \ D . Let u be the solutionof the following problem − ∆ u = f in Ω , u = g on Γ , u = 0 on ∂ Ω , (4.1)where g ( x, y ) = (4 − x )(4 − y ) on Γ and extended to Ω as ˜ g ( x, y ) = (4 − x )(4 − y ) cos(1 − x − y ) . Setting the solution equal to: u ( x, y ) = (4 − x )(4 − y ) (cid:0) χ D + exp(1 − x − y ) χ ¯ D (cid:1) , (4.2)where χ D i , i = 1 , f is chosenas: f = (cid:40) − ∆ u in Ω \ Γ , . All the computations have been performed employing a Voronoi dual mesh of a Delaunay trian-gulation (i.e., the dual mesh is obtained by connecting the barycentres of the triangles with straightlines).To validate the estimates (3.4) we consider in a separate way the influence of h and (cid:15) on the error.More precisely, we first explore the convergence with respect to h and then we study the convergencewith respect to (cid:15) . In both cases we consider a uniform discretization of the domain Ω so to have κ = δ = h . Convergence w.r.t. h . We set (cid:15) = 2 − (cid:28) h while we let h vary as h = 0 . , . , . , . . From Figure 3 we observe that the L -norm of the error decreases with order 1 while the errordecreases with order 1/2 in the H -norm. These rates of convergence are in agreement with (3.4). Remark 4.1.
If a local refinement of the diffuse interface region is performed in such a way that δ (cid:39) κ (cid:39) h (Figure 5), then first and second order of convergence are recovered for H and L norms,respectively (cf. [20, Section 6]). Convergence w.r.t. (cid:15) . We employ a fine mesh ( h = 0 . (cid:15) vary as: (cid:15) = 2 i , i = − , ..., − . (cid:15) (cf. (3.4))are obtained both in the L -norm (order 1) and in the H - norm (order 1 / (cid:15) becomes smaller than the chosen value of h , a plateau is observed as the(fixed) contribution from the discretization of the PDE (related to h ) dominates over the contributionfrom the introduction of the diffuse interface (related to (cid:15) ). e rr o r L2 error
Numerical orderOrder-1Order-1/2 0.01 0.02 0.03 0.05h10 H1 error
Numerical orderOrder-1Order-1/2
Figure 3: Error behaviour with respect to h (fixed (cid:15) = 2 − ): (left) L -norm error, (right) H -normerror. Dashed lines are theoretical convergence orders. e rr o r L2 error - h=0.006937
Numerical orderOrder-1Order-1/2 10 H1 error
Numerical orderOrder-1Order-1/2
Figure 4: Error behaviour with respect to (cid:15) (fixed h = 0 . L -norm error, (right) H -normerror. Dotted lines are theoretical convergence orders. In this paper we introduced a diffuse interface variant of a finite volume method, namely of the theso-called box method and obtained L and H error estimates highlighting the contributions from thediscretization parameter h associated to the polygonal computational mesh and the width (cid:15) of thediffuse interface. Despite the simplicity of the method, the present contribution seems to be novel inthe literature. Moreover, the present work may represent the first step towards the study of the diffuseinterface variant of more sophisticated finite volume schemes (possibly for more complex differentialproblems).This work opens fictitious boundary methods analysis to the box method and finite volume frame-work. Possible extensions of this research could be to being able to apply the plenty of penalizationmethods that are mostly thought for finite element implementations such as shifted boundary, Nitschepenalty, cut-fem or Brinkman penalization. 5 .200.100.050.01 h10 e rr o r L2 error
Numerical orderOrder-1Order-2 0.10h10 H1 error
Numerical orderOrder-1Order-2
Figure 5: On the left: example of a dual mesh with local mesh refinement around surrogate boundary.On the right: error behaviour with respect to h with local mesh refinement around the interface (fixed (cid:15) = 2 − ): (left) L -norm error, (right) H -norm error. Dashed lines are theoretical convergenceorders. The first author acknowledges the financial support of Fondazione Politecnico. The third authoracknowledges the financial support of PRIN research grant number 201744KLJL “
Virtual ElementMethods: Analysis and Applications ” funded by MIUR. The second and third authors acknowledgethe financial support of INdAM-GNCS.
References [1] Ivo Babuˇska. The finite element method with Lagrangian multipliers.
Numer. Math. , 20:179–192,1972/73.[2] Ivo Babuˇska. The finite element method with penalty.
Math. Comp. , 27:221–228, 1973.[3] Randolph E. Bank and Donald J. Rose. Some error estimates for the box method.
SIAM Journalon Numerical Analysis , 24(4):777–787, 1987.[4] Stefano Berrone, Andrea Bonito, Rob Stevenson, and Marco Verani. An optimal adaptive ficti-tious domain method.
Math. Comp. , 88(319):2101–2134, 2019.[5] Daniele Boffi and Lucia Gastaldi. A finite element approach for the immersed boundary method.volume 81, pages 491–501. 2003. In honour of Klaus-J¨urgen Bathe.[6] Christoph B¨orgers and Olof B. Widlund. On finite element domain imbedding methods.
SIAMJ. Numer. Anal. , 27(4):963–978, 1990.[7] Erik Burman and Peter Hansbo. Fictitious domain finite element methods using cut elements: I.A stabilized Lagrange multiplier method.
Comput. Methods Appl. Mech. Engrg. , 199(41-44):2680–2686, 2010.[8] Erik Burman and Peter Hansbo. Fictitious domain finite element methods using cut elements:II. A stabilized Nitsche method.
Appl. Numer. Math. , 62(4):328–341, 2012.[9] Philippe G. Ciarlet.
The finite element method for elliptic problems , volume 40 of
Classics inApplied Mathematics . Society for Industrial and Applied Mathematics (SIAM), Philadelphia,PA, 2002. Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58
Isogeometric analysis . John Wiley& Sons, Ltd., Chichester, 2009. Toward integration of CAD and FEA.611] J. Donea, S. Giuliani, and J.P. Halleux. An arbitrary lagrangian-eulerian finite element methodfor transient dynamic fluid-structure interactions.
Computer Methods in Applied Mechanics andEngineering , 33(1):689 – 723, 1982.[12] Richard E. Ewing, Tao Lin, and Yanping Lin. On the accuracy of the finite volume elementmethod based on piecewise linear polynomials.
SIAM J. Numer. Anal. , 39(6):1865–1888, 2002.[13] V. Girault and R. Glowinski. Error analysis of a fictitious domain method applied to a Dirichletproblem.
Japan J. Indust. Appl. Math. , 12(3):487–514, 1995.[14] Roland Glowinski, Tsorng-Whay Pan, and Jacques P´eriaux. A fictitious domain method forDirichlet problem and applications.
Comput. Methods Appl. Mech. Engrg. , 111(3-4):283–303,1994.[15] W. Hackbusch. On first and second order box schemes.
Computing , 41(4):277–296, 1989.[16] C. W. Hirt, A. A. Amsden, and J. L. Cook. An arbitrary Lagrangian-Eulerian computing methodfor all flow speeds [J. Comput. Phys. (1974), no. 3, 227–253]. volume 135, pages 198–216.1997. With an introduction by L. G. Margolin, Commemoration of the 30th anniversary { of J.Comput. Phys. } .[17] Thomas J. R. Hughes, Wing Kam Liu, and Thomas K. Zimmermann. Lagrangian-Eulerian finiteelement formulation for incompressible viscous flows. Comput. Methods Appl. Mech. Engrg. ,29(3):329–349, 1981.[18] X. Li, J. Lowengrub, A. R¨atz, and A. Voigt. Solving PDEs in complex geometries: a diffusedomain approach.
Commun. Math. Sci. , 7(1):81–107, 2009.[19] Charles S. Peskin. The immersed boundary method.
Acta Numer. , 11:479–517, 2002.[20] Matthias Schlottbom. Error analysis of a diffuse interface method for elliptic problems withDirichlet boundary conditions.
Appl. Numer. Math. , 109:109–122, 2016.[21] Jinchao Xu and Qingsong Zou. Analysis of linear and quadratic simplicial finite volume methodsfor elliptic equations.