Implicit Boundary Integral Method for Homogeneous Hele-Shaw Problem with multi-connected Domain

In this work, the implicit boundary integral method is implented to the homogeneous HeleShaw problem with a multi-connected domain. This method based on the solution of layer potential integral for the Laplace equation. The numerical technique is easy to implement, based on the idea of averaging the parameterization near the boundary and applying the Coarea formula. This technique changes the boundary integral into the Riemann integral that is easy to compute numerically. The difficulty in the computation of hypersingular integral occurs to compute the normal velocity of free boundary. The collocation technique is applied to eliminate the hypersingular part in the integral equation. Also, the numerical results are shown and its computation performance due to the appearance of a non-invertible matrix. In result, this method is a new technique for such a boundary problem with implementing level set method.


I. INTRODUCTION
A. Motivation I T is well known for the finite difference method (FDM) and finite element method (FEM) are usually used to solve a partial differential equation numerically.In this article, the boundary integral is considered to solve our problem, which is the equation for the evolution of the boundary equation.This method works based on the integral equation in the potential theory of elliptic problem.In Hale-Shaw problem, as seen in (8) free boundary problem appears in solving the elliptic equation.The function u is defined in the time depended domain t , where its boundary evolving caused by the gradient of u and some positive function g.Computationally, to solve the equation at time t, the domain at time t must be known.FDM is simply to be implemented, but computationally it is not efficient for moving boundary, since the domain mesh grid is not fixed in time.It happens as well in FEM; this method also depends on constructing the mesh grid for a given domain.If the given domain change in time, then in implementing FEM it is necessary to reconstruct the mesh grid in every time step.It causes to a higher cost of computation time.Boundary integral method (BIM), considers only the change near the boundary.Furthermore, the solution is constructed in the term boundary integral by finding a density defined on the boundary by fitting a given boundary condition.These ideas relate to surface potential integral can be read further in [4], [3], [1].In [5], they propose an implicit way to present the boundary integral equation using Coarea formula, such that the surface integral turn into Riemann integral equipped with its signeddistance representation.An important advantage of implementing this method that explains in [5] is this method is reliable for any arbitrary domain.This sustains a necessity for moving boundary problem with general velocity movement.Hele-Shaw problem is known as a popular model of incompressible liquid's flow in between two parallel plates called Hele-Shaw cell.Such a cell firstly was introduced by Henry Selby Hele-shaw as a practical model shown in his class to model a pressure driven flow of the fluid is injected into a shallow geometry from the center of two parallel plates when the surface of the fluid is bounded by another different viscous fluid or the air.The Hele-Shaw cell concept is also implemented to model dendritic form in crystalization, which has been studied following the growth pattern of Safman-Taylor in Hele-Shaw cell [7].Also, this model is used in the application of the surface evolution between viscous fluid and melting crystal, for instance, solidification or melting process.In three dimensions, this model is used as a model of the porous medium problem, for instance the flow of water absorption inside of the sand.In the industrial application, it is used in the plastic molding, importantly to determine the way out for the melting plastic when it is filling the mold, which is the boundary between the melting plastic and the air.For the detail explanation how this model works for plastic molding and petroleum extraction, see [6].

B. Hele Shaw Problem
In this work, we perform Hele-Shaw problem in two dimensions, and the type of fluid we are working on is Newtonian viscous fluid.The pressure driven flow is caused by a given source with constant pressure in the middle of the parallel plates.The pressure driven flow in such phenomena is following the Darcy's Law that is a flow rate through the porous medium is proportional to the minus of pressure gradient of viscous fluid.Accurately, it also depends on the viscosity of the fluid and permeability of the media, but we assume they are constant for simplicity, such that from conservation law of mass, it is obtained the equation for the density changing in time as follows, where ⇢ denotes the density over the unit volume, and u is a pressure.Since we assume the fluid is incompressible, then the pressure satisfies only the Laplace equation.The governing equation can be derived from the three dimensions Navier-stock equation by averaging the gap.Let 0 be a bounded open domain.We assume that in time t the fluid is filling the domain t R 2 , and the boundary t is a boundary between the liquid and the air.When the injected fluid evolves slow enough and the height h of the plates is also small, the averaged velocity q over of gap satisfying Since we assume the permeability and the viscosity to be constant, without losing a generality, the averaged of velocity proportional to the negative of pressure gradient, such that q = u.Therefore, the flow is almost stationary in between the gap of the plates.
where K is the area inside of t that is injected by the water continuously.There are two boundary conditions in the moving boundary t , the kinematic and the dynamic boundary condition.Since this Irma Palupi Implicit Boundary Integral... 94 work is neglecting the surface tension, the dynamic boundary condition is given to be, For the kinematic boundary condition, the fluid remains on the boundary in every time t with normal velocity V n , where n is the outer normal velocity.Since the pressure in ( t ) c is set to be zero (the water pressure in the air), u can be understood as a level set, such that Such a problem is called homogeneous Hele-Shaw problem with the media representation is given as a constant, that is homogeneous all over the surface of the plates.Furthermore, instead of putting the media as a constant, it modifies to be a positive continuous function depends on position and time, let us denote as a function g(x, t), where 1 g(x,t) can be understood as a depth of the hole must be filled by the fluid when the domain is evolving.Let K be a compact set with a smooth boundary satisfying 0 K, where 0 is an open domain in R n with also smooth boundary.
For given a positive continuous function g(x, t), closed set K, and an initial domain 0 K, the Hele-Shaw problem is reformulated as the problem of finding the pair solution of (u, t ) satisfying ( 8), with initial condition u(., 0) satisfying the Laplace equation in ( 0 \ K) with suitable boundary data as in (8).

II. LEVEL SET METHOD AND NUMERICAL ALGORITHM
Fast sweeping method is implemented in order to determine a level set function.Given a closed hypersurface , it is evolved by giving only the velocity in the normal direction.A boundary is defined through level set representation as a zero level set of level set function F , i.e for each point y , it satisfies F (y(t), t) = 0. Thus, the level set evolution is following equation (9).
where V n := n y • y (t) is denoted as a normal velocity.

A. Initialization Technique
In the level set representation, : {y : F (y) = 0}, and its normal vector n y is defined as Therefore, respectively for the interior and exterior closest point x the distance to the boundary can be computed exactly as finding d > 0, such that C, the distance of x is computed by using where n(x) = If it is given only a level set data, the interpolation to compute the distance of the closest points to the boundary may be needed, since the position of the boundary curve does not intersect the grid points.Figure .1 describes the idea how to interpolate the distance at node in 2-D rectangular mesh grid.This simple idea can be used in a case the curve is nice enough being in C 1 ( ).Given a level set data F (x i,j ) on the meshgrid, it assumes the function F is piecewise linear function in the domain, such that if the multiplication of F (x i,j ) and the level set in its neighbors points returns negative, then the zero level set {x|F (x) = 0} exists along the edge connecting x i,j to its neighbors.Then, the distance of x i,j to the boundary is computed as where k indicates the surrounding points on the hexagon in Figure . 1 with x i,j as a center.

B. Fast Sweeping Algorithm
The notion of Fast Sweeping algorithm based on Godunov upwind difference scheme to solve Eikonal equation.A detail description about this method can be found in [9].In the case of distance function, a viscosity solution d(x) 0 and it satisfies the following Eikonal equation.
Problem ( 12) is discretize using Godunov upwind difference scheme, such that partial derivative of d(x), for x = (x 1 , x 2 , . . ., x n ) R n is approximated to be Therefore, for n = 2, x the discrete form for (12) can be written as where, di = min{d i+1,j , d i 1,j } and dj = min{d i,j+1 , d i,j 1 }.By applying a periodic boundary condition, the one-sided nodes satisfy the following situation.
Gauss-Seidel iteration to sweep the whole distance value on the meshgrid.As describe in [9], it applied the computation for rectangular grid through four alternating directions as the following ordering, such that for every (i, j) do the computation as follow.

III. IMPLICIT BOUNDARY INTEGRAL METHODS
In level set method, the interface is described as a level set function in a zero level set.In a problem (9), the interface can be defined as {x |F (x) = 0}, where F is a level set function.One of level set function that commonly used is a signed distance function.A signed distance function can be used to describe implicitly the interface, and it has a gradient that does not vanish around .Moreover, it has a nice property that is | d| = 1.Distance function is a level set function that returns the shortest distance to the nearest point on the interface.In the interface which usually can be described as a closed curve (in 2D) or closed surface, it uses a Signed distance function to sign whereas the point in interior or exterior of domain surrounded by the interface.Signed distance function respect to boundary is defined as follow.
One of the important properties of signed distance function if the boundary sufficiently smooth is, the distance function is also smooth in some tubular neighborhood of , linear and has slope 1 along the normal direction to the interface, i.e | d| = 1.Therefore, by using d(x) as a level set function, the unit outer normal vector of can be obtain by the following formula.
According to the form of boundary integral, signed distance function is required to represent the boundary in the integral formula implicitly.The idea of the formulation is described by [5].The notion is based on averaging the parameterization of boundary integral equation by using delta Dirac function, then use the Coarea formula, such that the boundary integral changes the form into the Riemann integral defined in R n .The signed distance function lies in the new form of integral form as we compute the value.Regularization of delta Dirac function turnout the computation consider only the points on the tubular neighborhood.See Figure .2 for the illustration.From this state, the well known numerical integral technique can be used for computation.

A. Solution in the form of boundary integral
For simplicity, denote t = K t .The solution of ( 8) can be constructed as a boundary integral form.We consider Double-Layer Potential (DLP) integral to construct the solution of Laplace problem in t \ K, which is part of Hele-Shaw.Then, for x \ K the solution of ( 8) is: where is the fundamental solution of Laplace equation defined as the following function.
(x, y) =  The equation ( 18) is fundamental solution of Laplace's equation, and for fixed y R n , it is a harmonic function in R n \ {y}.And the second derivative of (18) satisfies the following.
The unknown density function defined on the boundary can be determined by fitting the boundary data, such that it satisfies the limit of integral DLP for u(x) when x approaching t .But the boundary integral in equation ( 17) is only defined in interior of domain, and to fit the density function, the limit of DLP when it is approaching the boundary is applied as in the Theorem.(1) from [4].Theorem 1: For U of class C 2 , the double-layer potential with continuous density can be continuously extended from U to Ū and from Ū c to U c with limiting values and the integral exists as an improper integral.Therefore, according to problem (8), the density can be computed by solving the folowing boundary integral equation for x t .
where f (x) is just a boundary condition written in (8), Irma Palupi Implicit Boundary Integral... 98

B. Changing the Integral form
Consider the two following integrals.
(y(s)) (x, y(s))dy(s) (y(s)) (x, y(s))dy(s) Suppose that the arch length of a closed curve R n is 1, and let y(s ) : s [0, 1] be a parameterization on boundary , then for the continuous function f of our surface integral can be calculated as the following.
We consider the integral of function that lies on as an integral of potential.
Therefore, after using the chain rule, averaging the parameter with delta Dirac function, and applying the Coarea formula, now we obtain another form of DLP boundary integral (21) for a solution of Laplace this method the boundary is described as a level set function and is also grid independent, this method has a benefit for such a free boundary problem as considered in this work.The improvement is also necessary to increase the accuracy and the efficiency in the computation technique, specially in handling the ill-conditioned matrix that exist during computation.

Figure 1 :
Figure 1: Ilustration how to initialize a distance function at node (i, j), in case the boundary intersects two points of the node's neighborhood edges.

Figure 2 :
Figure 2: figure Tubular neighborhood of

Table I :
Error comparison Fast Sweeping computation performance between wo types of Initialization.