The immersed interface method for axis-symmetric problems and application to the Hele–Shaw flow

Applied Mathematics and Computation 264 (2015) 179–197

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The immersed interface method for axis-symmetric problems and application to the Hele–Shaw flow Juan Ruiz Álvarez a,∗, Zhilin Li b a b

Department of Physics and Mathematics, Universidad de Alcalá (UAH), Spain Center for Research in Scientific Computation (CRSC) and Department of Mathematics, North Carolina State University (NCSU), USA

a r t i c l e

i n f o

MSC: 65D17 65M06 65N06 Keywords: Immersed interface method Axis-symmetric interface problem Discontinuous coefficients Finite difference method Level set method Hele–Shaw flow

a b s t r a c t Many physical application problems are axis-symmetric. Using axis-symmetric properties, many three dimensional problems can be solved efficiently using two dimensional axissymmetric coordinates. In this paper, the immersed interface method (IIM) in axis-symmetric coordinates is developed for elliptic interface problems that have a discontinuous coefficient, solution or flux. A staggered grid is used to overcome the pole singularity. Other challenges include deriving the jump relations in axis-symmetric coordinates, designing the numerical algorithm when the interface is close to the pole (r = 0); computing interface quantities such as the normal and tangential directions, surface derivatives, curvature information, etc. The numerical algorithm is based on a finite difference discretization and uniform grid in the axis-symmetric coordinates. The finite difference scheme is the standard one away from the interface but is modified at grid points near and on the interface. The method is shown to be second order accurate in the infinity norm. The developed new IIM is applied to the Hele–Shaw flow and compared with the results from the literature. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Since its invention, the immersed interface method (IIM) [1,2] has become an efficient and accurate numerical method for solving interface problems. The IIM uses simple structured meshes and standard numerical schemes away from an interface. The numerical schemes are modified only near or on the interface to take into account the jump conditions across the interface so that the accuracy will not be deteriorated. The IIM has been applied to quite a few challenging problems. However, until recently, most of the research done on the IIM has been based on Cartesian meshes and on polar coordinates [3]. The purpose of this paper is to develop a second order sharp interface method for axis-symmetric elliptic interface problems with uniform meshes in axis-symmetric coordinates. Many physical application problems are axis-symmetric. Using axissymmetric properties, many three dimensional problems can be solved efficiently using two dimensional axis-symmetric coordinates. There is indeed a rich literature on different methods for interface problems based on structured meshes. Nevertheless, there are hardly second order accurate sharp interface methods for axis-symmetric problems. In [4], an immersed interface method is claimed to be developed for an axis-symmetric Stokes flow with little detail and accuracy check. We are aware that the IIM is being investigated in [5] for axis-symmetric Stokes flows, which has not appeared in the literature yet.

∗

Corresponding author at: Matemáticas, Departamento de Física y Matemáticas, Universidad de Alcalá, 28871 Madrid, Spain. Tel.: +34918856413. E-mail addresses: [email protected] (J. Ruiz Álvarez), [email protected] (Z. Li).

http://dx.doi.org/10.1016/j.amc.2015.03.131 0096-3003/© 2015 Elsevier Inc. All rights reserved.

180

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

Fig. 1. A diagram of the elliptic interface problem in axis-symmetric coordinates, where  = + − , with an interface  . The coefficient β may have a finite jump across the interface.

We will start with the elliptic interface problem

∇ · (β∇ u) − σ u = f, (r, z) ∈  \ 

(1)

in axis-symmetric coordinates. In the solution domain, we assume that there is an axis-symmetric interface   C2 across which the coefficients β (r, z) > 0 and σ (r, z), the source term f, the solution u, and the flux β un can have a finite jump, where n is the normal direction pointing outward, and un = u · n, see Fig. 1 for an illustration. To make the problem well-posed, often the jump conditions in the solution and the flux are known a priori

[u] = w,

[β un ] = v,

(2)

C2

C1

and v  are two functions that have co-dimension one compared with that of the solution u. where we assume that w  For the special case when w = 0, the above interface problem is equivalent to Peskin’s model

∇ · (β∇ u) − σ u = f (r, z) +





  vδ x(s) − X(s) ds,

(r, z) ∈ 

(3)

where s is some parametrization of the interface  . Note that the above PDE is defined on the entire domain including the interface  . There are a number of challenges in developing accurate immersed interface methods for axis-symmetric elliptic interface problems. First of all, it is not clear how to obtain the jump relations needed for the immersed interface method in axis-symmetric coordinates. Secondly, it is not clear how to deal with the pole singularity at the interface. Furthermore, it is not trivial to get the interface information needed for the IIM such as the normal and tangential directions, the curvature information, etc. These issues will be addressed in detail in this paper. In our numerical method, we use a staggered mesh in the r direction and the reflected boundary condition at the interface when it touches the pole. This is justified since it is reasonable to assume that the solution is piecewise smooth. We also show that the interface relations are almost the same as for the Cartesian mesh in two dimensions. The interface will be represented as the zero level set of a Lipschitz continuous function ϕ (r, z). We also applied our new IIM in axis-symmetric coordinates to a moving interface problem in tracking the location of the boundary in a Hele–Shaw cell. The application is challenging because it is unstable for some normal modes. A level set method in the axis-symmetric coordinates also needs to be developed including the re-initialization step, updating the level set function in a computational tube. Our method performs well against the results in the literature. Particularly, we have checked the results against the adaptive boundary integral method. The rest of the paper is organized as follows: In Section 2 we derive in detail the interface relations for elliptic interface problems and the numerical method in the axis-symmetrical coordinates. In Section 3, we show some numerical results and convergence analysis for some problems for which we know the analytic solutions. In Section 4, the newly developed IIM is

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

181

applied to a moving interface problem of the axis-symmetric Hele–Shaw flow. Some interesting numerical results are obtained. In Section 4, some conclusions and acknowledgements are given. 2. The numerical method Now we discuss the finite-difference method in axis-symmetrical cylindrical coordinates for a domain r1  r  r2 , z1  z  z2 where the parameters of Eq. (1) can have a finite jump across the interface, see Fig. 1 for an illustration. The goal is to numerically obtain an accurate solution u when the domain is divided through the interface into two different subdomains. The different subdomains will have different parameters. The elliptic equation shown in (1) can be written in axis-symmetrical cylindrical coordinates as:

1 ∂ r ∂r





r

∂β u ∂ 2β u + − σ u = f (r, z), (r, z) ∈  \  , ∂r ∂ z2

(4)

with a boundary condition at . The right hand side f is the source term. As stated before, the natural jump conditions across the interface  are

[u] = w(s),

[β un ] = v(s),

(5)

where s is a parameter used to represent the interface  (r(s), z(s)). The jump conditions mean that the solution and the flux can be discontinuous across the interface. Note that even if the flux is continuous across the interface  , the normal derivative may have a jump due to the discontinuity in the coefficient β . In this paper, we assume that   C2 . Also, for our practical application, we assume that the coefficient β is a piecewise constant,



β(x) =

β + if x ∈ + , β − x ∈ − ,

(6)

where x = (r, z). Even so, we provide the interface relations for more general cases, where β can be a function of r and z. We assume that the domain  is a regular one in axis-symmetric cylindrical coordinates, that is, : { (r, z), 0  r  R1 , Z0  z < Z1 } within which there is an arbitrary smooth interface. Note that our discussion is also valid and simpler when r1  r  R1 for r1 > 0. The interface is represented using the zero level set of a Lipschitz function ϕ = ϕ (r, z). To obtain an approximate solution of the elliptic interface problem, we set a uniform grid in cylindrical coordinates:

ri = r1 + i r, zj = z1 + j z,

M+1 , 2 j = 0, 1, · · · , N,

i = 0, 1, · · · ,

(7) (8)

for two given parameters M and N, with M being an odd number, where r = 2R1 /M and z = (Z1 − Z0 )/N. Thus, in our computations, r1 = h/2, z1 = h/2 and h is the step size in both directions. As shown in (4), the partial differential equation in axis-symmetric cylindrical coordinates is:

1 ∂ r ∂r



βr







∂u ∂ ∂u + − σ u = f (r, z). β ∂r ∂z ∂z

(9)

We use a finite difference method to discretize the partial differential equation in cylindrical coordinates. At a regular grid point (the grid points in the 5-point central finite difference stencil are all from one side of the interface  ) the discretized equation is



β

ui,j−1 − 2ui,j + ui,j+1 1 ri− 12 ui−1,j − (ri− 12 + ri+ 12 )ui,j + ri+ 12 ui+1,j + ri ( r)2 ( z)2



− σ ui,j = f (ri , zj ),

(10)

where ui, j is the finite difference approximation to the solution u(ri , zj ). As mentioned before, in order to deal with the pole singularity, we use a staggered grid. That means that the r coordinate of every grid point takes the value

ri = (i − 1/2) r,

r =

2R1 , M

i = 1, 2, · · · ,

M+1 . 2

(11)

Note that r1 = r/2 = h/2 and r M+1 = R1 , with M being an odd number. This configuration makes it easy to deal with the boundary 2

condition on the left side of the boundary. As shown in Fig. 2, if we set a ghost point at r = − r/2, then Eq. (10) for the regular points at r = r/2 takes the form

β



u1,j−1 − 2u1,j + u1,j+1 u−1,j − 2u1,j + u2,j 1 u2,j − u−1,j + + 0.5 r 2 r ( r)2 ( z)2



− σ u1,j = f (r1 , zj ),

(12)

where u−1, j is the value of the solution at r = − r/2, if we use the centered 5-point stencil. Simplifying, the result in (12), it turns out that the ghost point disappears from the equation



2

u1,j−1 − 2u1,j + u1,j+1 u2,j − u1,j + ( r)2 ( z)2



− σ u1,j = f (r1 , zj ).

(13)

182

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

Fig. 2. Example of an irregular point close to the left boundary of the domain and the points chosen for the stencil. The first column is composed by ghost points. .

Fig. 3. Example of an irregular point and the points chosen for the stencil.

Obviously, the finite difference scheme above is inaccurate at irregular grid points (i.e. the grid points in the centered 5-point stencil are from both sides of the interface), see Fig. 2 for a typical case. The first column of points depicted in Fig. 2 is composed by ghost points. As shown in (13) these points do not appear in the solution. The only exception is the point at (−1, j) in Fig. 2, where we set the reflected Neumann boundary condition ∂∂ur |r=0 = 0. This is justified since we assume that the solution is

piecewise smooth except at the interface. Otherwise there would be a discontinuity of ∂∂ur at r = 0. The reason is that the point at (1, j) is irregular, so we will need six points to find out the coefficients of the finite difference scheme. One of these points will be the one at (−1, j). We will discuss the IIM process used to construct the finite difference scheme at irregular grid points in the next subsection. 2.1. The finite difference scheme at irregular grid points The derivation of the finite difference scheme at an irregular point (ri , zj ) (see Fig. 3 for an illustration) for a discontinuous coefficient β using the IIM includes the following main steps: • Choosing a finite difference stencil. In our application, we choose six grid points that include the standard centered five-point stencil, plus one from one of the four corners, (ri ± 1 , zj ± 1 ). In particular, the closest one to the projection of the central point of the stencil over the interface. In order not to lose accuracy, we will choose a point that also belongs to the side of the interface where the central point falls. The finite difference scheme can then be written as

γ1 Ui−1,j + γ2 Ui,j + γ3 Ui+1,j + γ4 Ui,j+1 + γ5 Ui,j−1 + γ6 Ui+l,j+p − σij Ui,j = fij + Cij ,

(14)

where γ k are the finite difference coefficients that are going to be found using an undetermined coefficient method, l and p are either 1 or −1, and Cij is a correction term. The local truncation error at the grid point (ri , zj ) is

Tij = γ1 u(ri−1 , zj ) + γ2 u(ri , zj ) + γ3 u(ri+1 , zj ) + γ4 u(ri , zj+1 ) + γ5 u(ri , zj−1 ) + γ6 u(ri+l , zj+p ) − σij u(ri , zj ) − fij − Cij . (15) • Choosing a point on the interface (r∗ , z∗ ) that is close to the irregular grid point (ri , zj ). Often the point is chosen as an approximation of the orthogonal projection of (ri , zj ) on the interface. In Section 2.3 we give a brief explanation on how to obtain the orthogonal projection. In the neighborhood of (r∗ , z∗ ), the interface can be written as ξ = ξ (η) with ξ ∗ = ξ (η∗ ). • Setting up the local coordinates system at (r∗ , z∗ ). To do so, we assume that the interface can be expressed in terms of the local coordinates as  : { (ξ , η), ξ = ξ (η)}, that is a closed curve. As mentioned before, this curve is approximated using the zero level set ξ = χ (η) of the function ϕ = ϕ (ξ , η). Thus, the interface can be represented as the zero level set ϕ (ξ , η ) = χ (η ) − ξ .

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

183

• Expanding u(ri + l r, zj + p z) at (r∗ , z∗ ) from each side of the interface, where l and p take the value 0, −1, or 1 ∗ ± u(ri + l r, zj + p z) = u± (r∗ , z∗ ) + (ri + l r − r∗ ) u± r + (zj + p z − z )uz +

+

(zj + p z − z∗ )2 2

(ri + l r − r∗ )2 2

u± rr

∗ ∗ ± 2 u± zz + (zj + p z − z )(ri + l r − r )urz + O(h ),

where all the partial derivatives are evaluated at (r∗ , z∗ ) and r = z = h for simplicity. • Using the interface relations to represent all the quantities from one particular side in terms of the other side, for example:

u+ = u− + w  1 − − u+ = β u + v ξ ξ +

β

− u+ η = uη + wη ··· ··· ···

···

All the interface relations are derived in the next section. • Collecting terms in the expression for the local truncation error and setting-up a system of equations for the finite difference coefficients from the PDE at (r∗ , z∗ ) from the chosen side, say from the “−” side: − − − Tij = L1 (γ1 , . . . , γ6 )u− + L2 (γ1 , . . . , γ6 )u− ξ + L3 (γ1 , . . . , γ6 )uη + L4 (γ1 , . . . , γ6 )uξ ξ + L5 (γ1 , . . . , γ6 )uξ η   dw d2 w dv − + L6 (γ1 , . . . , γ6 )u− , , v, , ηη − fij − σij u − Cij + L7 γ1 , . . . , γ6 , w, dη dη 2 dη

(16)

where all Ll , l = 1, 2, , 7, are linear functions of their arguments. The system of equations for the coefficients is then:

L1





γ1 , . . . , γ6 = 0

cos θ − β + βξ− r sin θ − L3 (γ1 , . . . , γ6 ) = − β + βη− r L4 (γ1 , . . . , γ6 ) = β − L2





γ1 , . . . , γ6 =

(17) (18)

L5 (γ1 , . . . , γ6 ) = 0 L6 (γ1 , . . . , γ6 ) = β − . • Determining the correction term Cij to minimize the local truncation error in the magnitude:



Cij = L7

γ1 , . . . , γ6 , w,

dw d2 w dv , , v, dη dη 2 dη



.

Note that the process above is necessary for the IIM method at an irregular grid point. The interface relations needed to represent the quantities from “+” side in terms of “−” side will be derived in the next subsection. Also note that we will run into another problem if the irregular grid point is near the pole r = 0) as we can see from Fig. 1. As mentioned at the end of Section 2, in these cases we use a ghost point for every irregular point that we find on the boundary and set a reflected Neumann boundary condition ∂∂ur |r=0 = 0 at those points. 2.2. Deriving the interface relations in axis-symmetric cylindrical coordinates As we mentioned in the previous subsection, given the jump conditions, [u] = w and [β un ] = v, we need to express the quantities from one side of the interface in terms of the other side up to second order. Referring to Figs. 3 and 1, the interface  divides the domain into two parts, which we simply denote as “+” and “−” sides. Assume that the interface  at a particular point (r∗ , z∗ ) can be represented in terms of the local coordinates as ξ = χ (η), where η is the tangential direction τ and ξ is the normal direction n shown in Fig. 1. As mentioned in [2], we can express r and z in terms of the local coordinates as

 



ξ sin θ cos θ = − sin θ cos θ η



r − r∗ z − z∗



,

(19)

where θ is the angle between the x-axis and the normal direction. So we can define the local transformation matrix as



cos θ − sin θ

A=

 sin θ . cos θ

(20)

From (19) it can be easily obtained that the first order derivatives in the local coordinate system can be expressed as



uξ uη





=

cos θ − sin θ

sin θ cos θ

  ur , uz

(21)

184

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

and that the second order derivatives can be expressed as

⎞ ⎛ 2 cos θ uξ ξ 2 ⎝uηη ⎠ = ⎜ ⎝sin θ uξ η − sin θ cos θ ⎛

sin θ cos2 θ sin θ cos θ 2

⎞⎛ ⎞ 2 sin θ cos θ urr ⎟ −2 sin θ cos θ ⎠ ⎝uzz ⎠ . 2 urz cos2 θ − sin θ

(22)

It is reasonable to assume that u(x) is a piecewise C2 function in the neighborhood of  , excluding  . Below is the sketch of the proof of how to obtain the following interface relations.

u+ = u− + w

(23)

β + u+ξ = β − u−ξ + v

(24)

− u+ η = uη + wη

(25)

− + − u+ ηη = uηη + (uξ − uξ )χηη + wηη   β− β + u+ξ η = −βη+ + + βη− u−ξ + χηη (−β − + β + )u+η + β − u−ηξ + β + χηη wη −

β

(26)

βη+ v + vη . β+

(27)

Eqs. (23) and (24) represent the original jump conditions on the interface. If we differentiate (23) in terms of η, we obtain + − − u+ η + uξ χη (0) = uη + uξ χη (0) + wη .

(28)

Knowing that in a neighborhood of the interface ξ = χ (η) and χ (0) = χ η (0) = 0, we obtain Eq. (25). If we differentiate (25) in terms of η we obtain + + + + − − − − − u+ ηη + uηξ χη + (uξ η + uξ ξ χη )χη + uξ χηη = uηη + uηξ χη + (uξ η + uξ ξ χη )χη + uξ χηη + wηη .

(29)

Using again the fact that χ (0) = 0 and χ η (0) = 0 at the origin of the local coordinates system, we obtain Eq. (26). Before differentiating (24) in terms of η, we express the normal unit vector of the interface as

(1, −χη ) n=  , 1 + χη2

(30)

using the fact that the zero isocontour of the level-set function can be expressed as ϕ (ξ , η) = ξ − χ (η). Thus, we can write Eq. (24) as

⎡

⎤

β(uξ − uη χη ) ⎦ ⎣  1 + χη2

= v(η).

(31)

Differentiating the identity above in terms of η, we get

       χ χ βη + βξ χη uξ − uη χη + β uξ η + uξ ξ χη − uηη + uηξ χη χη + uη χηη = vη 1 + χη2 + v  η ηη . 1 + χη2



(32)

At (ξ , η) = (0, 0) we obtain Eq. (27). If the parameter β is constant, we can get rid of the terms multiplied by β η , β ξ . Using the differential equation (9), we can find the relation for u+ ξξ



1 ∂ r ∂r



βr





∂u ∂ ∂u + β ∂r ∂z ∂z



− [σ u] = [f (r, z)] .

(33)

Expanding the expression above we obtain



β ur r∗

 + βr ur + β urr + βz uz + β uzz − [σ u] = [f (r, z)] ,

(34)

where r∗ is the r coordinate of the projection on the interface. We can simplify this expression using the local coordinates system, where the projection over the interface is the origin of the system (ξ ∗ , η∗ ) = (0, 0). We can express the previous relation in terms of the local coordinates, knowing that



uξ ξ uηξ

uξ η uηη





=A

urr uzr

 urz AT , uzz

(35)

where A is given by (20). The expression (34) in the local coordinates is



β cos θ r

    −β sin θ + βξ uξ + + βη uη + β uξ ξ + β uηη − σ u = [f ]. r

(36)

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

185

If β is piecewise constant, the above relation is reduced to

 

β

 1 uξ cos θ − uη sin θ + uξ ξ + uηη r



= [f ].

The solution of the system of Eqs. ([23–27,36]) provides us with the expressions for the partial derivatives from both sides of the interface, i.e. the interface relations:

u+ = u− + w u+ ξ =

1

β+

(37)

(β − u−ξ + v)

(38)

− u+ η = uη + wη

u+ ηη = +

uξ ξ =

χηη

(39)

χηη v β+ − β− − uξ + u− ηη + wηη − β+ β+ 

1

β +r



rσ [u] −

β rχηη − rχηη β + +

−

(40)

rβξ+ β −

β+

− rβξ

−

u− ξ

 −  + − − − − (rβη+ + β − sin θ − β + sin θ − rβη− )u− η + r β uξ ξ − r β − β uηη  rβξ+ v + − v cos θ − rχηη v − β + wη sin θ + β + rwηη + + r β w − [f ]r η η + u+ ξη =



1



βη− −

β+

β

(41)



βη+ β − − βη+ v   − − + uξ + χηη β + − β − u− + vη . η + β uξ η + β χηη wη − + β β+

(42)

When σ = 0 and β is piecewise constant, we obtain the following interface relations:

u+ = u− + w u+ ξ =

1

β+

(β − u−ξ + v)

− u+ η = uη + wη

χηη v β+ − β− − uξ + u− ηη + wηη − β+ β+

u+ ηη =

χηη

u+ ξξ =

(−(β + rχηη − rχηη β − )u−ξ − sin θ (β − − β + )u−η + rβ − u−ξ ξ − β +r + + − r(β + − β − )u− ηη (v cos θ − r χηη v − β wη sin θ + β rwηη + wη − [f ]r ))

u+ ξη =

(43) (44) (45) (46)

1

1

β+

(χηη (β + − β − ) u−η + β − u−ξ η + β + χηη wη + vη ).

(47) (48)

Thus, we have obtained all the six interface relations needed for deriving a second order accurate finite difference scheme using the IIM. The scheme is exact if the solution is a piecewise quadratic function and the interface and its derivatives can be computed exactly. The derivation of the “−” quantities in terms of the “+” side is similar. 2.3. Obtaining the normal direction and the distance to the interface using the level set function We need to find the orthogonal projection of irregular grid points on the interface. As we solve problems in cylindrical coordinates, it is natural to obtain the normal direction in cylindrical coordinates. The cylindrical coordinates system can be expressed as

x = r cos θ ,

(49)

y = r sin θ ,

(50)

z = z,

(51)

with r  0 and 0  θ < 2π . As stated in [2], if x = (r, z) is an irregular point, the orthogonal projection x∗ on the interface can be written as,

x∗ = x + α p,

186

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

where p is the normal direction to the interface and α is a scalar. If the normal direction is known, then we can use the following expression to find α ,

1 2

ϕ(x) + (∇ϕ(x) · p)α + (pT He(ϕ(x))p)α 2 = 0.

(52)

This expression is the Taylor expansion of the level set function at the point x (i.e. the central point of the stencil in our case) around the projection x∗ = (r∗ , z∗ ). Thus, He(ϕ (x)) is the Hessian matrix of the level set function ϕ(x),

He(ϕ(x)) =





ϕrr ϕrz . ϕzz ϕzr

(53)

Solving the quadratic Eq. (52), we obtain a third order approximation to the distance from the central point of the stencil to the interface. The normal direction is given by the gradient of the level set function ϕ (x) = (ϕ r , ϕ z ). Once we have the projection x∗ , it is straightforward to set up a local coordinate system centered at x∗ . The unit normal direction will be given by,

ξ=

∇ϕ (ϕr , ϕr ) =  . |∇ϕ| ϕr2 + ϕz2

 z , −ϕr ) . Then, ξ and η form a local coordinate system in the normal and tangential The tangent direction is taken as η = (ϕ

ϕr2 +ϕz2

directions. 2.4. Obtaining the surface derivatives in 2D In order to obtain the transformations proposed in (37–42) it is necessary to obtain the surface derivatives χ η (which is zero at the projection) and χ ηη of the interface at the projection x∗ . These quantities can be easily obtained using the level set function. We know that,

ϕη + ϕξ χη = 0

(54)

ϕηη + ϕηξ χη + (ϕξ η + ϕξ ξ χη )χη + ϕξ χηη = 0.

(55)

Thus, knowing that χ (0, 0) = 0 and χ η (0, 0) = 0, we obtain

χη = −

ϕη ϕηη , χηη = − . ϕξ ϕξ

(56)

The quantities ϕ ξ , ϕ η and ϕ ηη can be obtained from (21) and (22) in terms of the original axis-symmetric coordinates at the grid points using second order central finite differences. As the projection typically does not fall at a grid point, we can use the bilinear interpolation to obtain these quantities at the projection and then substitute at expression (56) to obtain the surface derivatives. Another option is to work in the local coordinates and to consider the interface as a plane curve in the plane zr given in the implicit form ϕ (ξ (η), η) = 0. In this case, in the local coordinates system, we can obtain the curvature as [6],

κ =−

ϕηη ϕξ2 − 2ϕη ϕξ ϕηξ + ϕξ ξ ϕη2 (ϕη2 + ϕξ2 ) 2 3

.

As before, we can obtain the curvature at grid points using central differences and then use the bilinear interpolation to compute the curvature at the projection. 2.5. Computation of the derivatives of interface quantities If the jump conditions [u] = w and [β un ] = v are different from zero, then we will need the first order derivatives wη , vη and the second order derivative vηη in order to evaluate the transformation proposed in (37–42) and (43–48). These quantities can be obtained through the least squares interpolation, proposed in [7] and [2]. We know that v and w are functions defined only on the interface. Therefore, we suppose that we know their values at every projection {X∗p }. If we can obtain the values of v and w at several points near the projection X∗ , it is possible to obtain a nearly second order approximation of the derivatives of the interface quantities at the projection X∗ . We will explain the process using the least squares interpolation. In a neighbourhood of X∗ the general function g(X), defined only on the interface, can be written as g(η) in terms of the local coordinates (19) centered at X∗ . The least squares interpolation of gη and gηη can be expressed as,

gη (X∗ ) =



|X∗ −X∗p |≤R , p≤Np

αp gp

(57)

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197



gηη (X∗ ) =

λp gp ,

187

(58)

|X∗ −X∗p |≤R , p≤Np

where gp is the value of the function at each projection { X∗p }. Rε is a parameter that determines how many points are being used in order to obtain the interpolation. We have chosen this parameter so that Np = 12, i.e. we choose the 12 closest projections of irregular grid points over the interface. In order to improve the accuracy of the estimation of the derivatives gη , gηη , we can use projections only from one side of the interface. Now we are ready to explain how to obtain the values of α p for gη and λp for gηη . We will use the Taylor expansion and the singular value decomposition (SVD) in order to solve an undetermined system of equations. Using the local coordinates (19) centered at X∗ , we denote the η coordinate of an orthogonal projection point X∗ p as ηp . If we use now the Taylor expansion at X∗ , we obtain

 |X∗ −X∗p |≤R

|X∗ −X∗p |≤R

= where get



αp gp =



1 2

∗ αp g∗ + gη∗ (X∗p − X∗ ) + gηη (X∗p − X∗ )2 + · · ·

∗ ( )g∗ + ( )gη∗ + ( )gηη + ···

g∗ , g∗ , g∗ η



ηη are the function value and its derivatives at the point

αp = 0,

p





αp ηp = 1,

p

for gη∗ , and

λp = 0,

p







(59) X∗ .

Using the method of undetermined coefficients, we

αp ηp2 = 0,

(60)

λp ηp2 = 1,

(61)

p

λp ηp = 0,

p

 p

∗ . In these expressions, η are the coordinates of the projections X∗ in the parametric form ξ = χ (η ). The systems of for gηη p p Eqs. (60) and (61) are underdetermined (mind that in our case, p = 12) so the SVD should be used in order to obtain a pseudoinverse. The SVD provides the solution with the least 2-norm among all possible. It is important to realize that in order to obtain the different derivatives, the coefficient matrix is always the same, i.e. the left hand side of Eqs. (60, 61) is the same. Thus, we need to compute the SVD only once. With the least squares interpolation it is possible to obtain a nearly second order accurate gradient of a computed solution that is second order accurate. An alternative way of computing the derivatives of interface quantities can be found in [2].

3. Numerical examples and validations for elliptic interface problems with a fixed interface In this section we present some examples for elliptic interface problems with a fixed interface in axis-symmetric coordinates. Example 1 In the first experiment we set the coefficient β as a piecewise constant. We try to solve the Poisson problem  · (β u) = f in an axis-symmetric cylindrical domain. The jumps [u] in the solution, in the flux [β un ], and in the source term [f] are determined from the exact solution

u(r, z) =

 (r2 ) e · sin(z) if r 3 + z3 if

ρ ≤ 0.4 + 0.04 cos(6θ ), ρ > 0.4 + 0.04 cos(6θ ),

(62)

 were ρ = r2 + z2 and θ = arctan(z/r). We can see that the interface is a rather complex function: a flower of six petals (three petals in the interval θ  [ − π /2, π /2]). The domain used is 0  r  0.7, −0.7  z  0.7. The interface and domain described are precisely the one depicted in Fig. 1. For this example, both the solution and the normal derivatives are discontinuous across the interface. Table 1 shows the results of a grid refinement analysis. We list errors for L1 , L2 , and L norms. The approximate order of convergence is for the L norm. Note that, the jump ratio β − /β + is quite large/small in our test, which shows that our developed algorithm can handle large or small jump ratios. Example 2 In the second experiment we try to use the same geometry of the interface, but with the exact solution

u(r, z) =



r 2 + z2 if r2 + 3z + r + 1 if

ρ ≤ 0.4 + 0.04 cos(6θ ), ρ > 0.4 + 0.04 cos(6θ ).

(63)

The grid refinement analysis for this experiment is presented in Table 2 and Fig. 4. In both cases we can observe second order accuracy.

188

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197 Table 1 A grid refinement analysis for the solution presented in (62) with a different jump ratio in β − and β + . The numerical results are second order accurate in the L norm.

β + = 1000, β − = 1

n

32 64 128 256 512 1024

β + = 1, β − = 1000

||EN ||

||EN ||1

||EN ||2

Order

||EN ||

||EN ||1

||EN ||2

Order

1.1489e−02 1.7932e−03 1.8823e−04 3.3638e−05 8.4729e−06 1.6272e−06

1.6670e−03 2.3538e−04 2.5036e−05 5.3614e−06 1.4438e−06 2.7086e−07

6.4833e−06 1.1730e−07 1.4669e−09 6.2127e−11 4.5593e−12 1.5393e−13

– 2.6192 3.2613 2.4241 1.9835 2.1586

2.1896e−01 9.5442e−02 2.8332e−03 3.0183e−04 6.8734e−05 5.4514e−06

9.3765e−02 4.5056e−02 1.2186e−03 1.0686e−04 3.0342e−05 1.7628e−06

1.5583e−02 3.3711e−03 2.3942e−06 1.9010e−08 1.4679e−09 5.3293e−12

– 1.1709 5.0169 3.2124 2.1286 3.6512

Table 2 A grid refinement analysis for the solution presented in (63) with a different jump ratio in β − and β + . The numerical results are second order accurate in the L norm. n

32 64 128 256 512 1024

β + = 1000, β − = 1

β + = 1, β − = 1000

||EN ||

Order

||EN ||

Order

2.8917e−02 5.7908e−03 2.0948e−03 3.1037e−04 6.3054e−05 1.0516e−05

– 2.2677 1.4504 2.7393 2.2928 2.5803

2.317 1.1475 2.8216e−02 8.9452e−03 2.2269e−03 7.5260e−04

– 0.9908 5.2856 1.6480 2.0004 1.5629

−3

1 y=2.4555x +8.4794 0

−5

−1

−6

−2

−7

−3

log(|e∞|)

−4

∞

log(|e |)

y=2.2509x +3.5749

−8

−4

−9

−5

−10

−6

−11

−7

−12 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

−8 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

Fig. 4. Linear regression plot for the solution presented in (63) when β + = 1000 and β − = 1 (left) and when β + = 1 and β − = 1000 right. In order to obtain this plot, we have used the data shown in Table 2. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

In the previous experiments we can see that when studying the accuracy against the exact solutions, the order of convergence for β + = 1 and β − = 1000 varies if we compare it to the case β − =1 and β + = 1000. The answer to this fact is that the error depends on quite a few factors: the solution and the coefficient on each side, the curvature, the relative positions between the interface and the underlying grid . . . . Even so, the error behaves like C(κ , β + , β − ,  , N)h2 log (1/h) where C is of order one but not a constant. Note that the errors with different jump ratios are comparable when the grid is fine enough. 4. The level-set method for the Hele–Shaw flow in axis-symmetric coordinates In this section we will study the behavior of the Hele–Shaw flow in which the interface is moving. This flow is typically found between two closely spaced parallel glass plates [8]. Usually, the fluids used are immiscible and one has higher viscosity than

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

189

the other. We will put a sink in the viscous fluid to model the situation in which a less viscous fluid drives a more viscous fluid. The jump in the pressure across the interface that separates the two fluids is balanced by the surface tension. In this problem, the velocity field is proportional to the gradient of the pressure p. The governing equations are

u = −β∇ p,

(64)

∇ · u = φ ,

(65)

2

b with β = 12 μ , where b is the distance between the two plates and μ is the viscosity of the fluid, which is very different inside and outside the interface that separates the fluids. For the interested reader, there is a relevant literature about the physical behavior [9–13] and the numerical simulation of the Hele–Shaw flow [11,14–19]. In particular, this problem has been widely used as benchmark problem for numerical algorithms to compute unstable fronts. The problem can be written as,

∇ · (β∇ p) = −φ .

(66)

The source term is the result of the absorption of the most viscous fluid out of the Hele–Shaw cell in order to force the less viscous fluid to drive the more viscous fluid. We assume that the interface is far away from the boundary, in order to obtain the boundary condition for the pressure. The flow at the boundary agrees with the radial outflow produced by the source term in a uniform fluid. In our case, we have used,

p0 p(r, z) = − √ 2 r + z2 at the boundary, where p0 is an arbitrary constant. In our experiments we have used a delta function as source term. It is also possible to smear out the source term as shown in [20] in order to speed out the simulations. As can be seen, the problem can be transformed into a Poisson problem. Once we have obtained the pressure through the IIM, we can obtain the velocity using Eq. (64). Once the velocity has been obtained, it can be used to compute the normal velocity. This velocity can be plugged in the Hamilton–Jacobi equation to update the position of the interface. 4.1. Validation of the numerical method using a benchmark problem for the Hele–Shaw flow The Hele–Shaw flow has served as a benchmark problem for numerical algorithms, as is well known that it is intrinsically unstable. In the next example, we check the accuracy of our method for an exact solution of the Hele–Shaw flow. We will test the result for the pressure,

⎧   V0 2r3 3r2 ⎪ ⎪ ⎪ − + C1 ⎪ ⎪ β − 3α 2 2α ⎪ ⎪ ⎪ ⎨ V α p(ρ) = − 0− log(ρ) + C0 ⎪ β ⎪ ⎪ ⎪ ⎪ ⎪ V ⎪ 0α ⎪ ⎩− + log(ρ)

β

if 0 < ρ ≤ α , if

α < ρ < r ,

(67)

if r ≤ ρ ,

 with ρ = r2 + z2 and r = 0.5 the radius of the initial interface. In our experiments we have used V0 = 1, β + = 1, β − = 100 and α = 0.1. Note that the pressure is continuous across the circle r = α . C0 is chosen as

C0 =

γ r

+ V0 α log(r )



1

β−

−

1

β+



,

such that [p]r = γ κ = −γ /r . Where γ is the surface tension and κ is the curvature of the interface, that arise from the Laplace–Young condition. The pressure is continuous across the inner boundary r = α . The constant C1 ,

C1 = C0 −

V0 α log(α)

β−

+

5V0 α , 6β −

is chosen to satisfy the condition [β pn ]r = α = 0. The pressure p satisfies the following Poisson equation

∇ · (β∇ p) =

  ⎧ −8ρ + 9α ⎪ ⎪ −V if 0 < ρ ≤ α , 0 ⎨ α2 ⎪ ⎪ ⎩−V0

α ρ2

if

(68)

α < ρ.

The grid refinement analysis for this problem is presented in Table 3 and in Fig. 5. In both cases, the results show second order accuracy. This example also validates our simulation for the Hele–Shaw flow.

190

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197 Table 3 A grid refinement analysis for the solution presented in (67). The order in the table is for the L norm. The numerical results are second order accurate for the L norm. n

32 64 128 256 512 1024 2048

β + = 1, β − = 100 ||EN ||

Order

2.7729e−03 3.9168e−04 1.0750e−05 7.5266e−06 2.4778e−06 9.1752e−07 2.4330e−07

– 2.7599 5.1288 0.5113 1.5984 1.4312 1.9137

−4 y=2.1311x +−0.96595 −6

∞

log(|e |)

−8

−10

−12

−14

−16 −7

−6.5

−6

−5.5

−5

−4.5 log(h)

−4

−3.5

−3

−2.5

Fig. 5. Linear regression plot for the solution presented in (67) when β + = 1 and β − = 100. In order to obtain this plot, we have used the data shown in Table 3. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

4.2. Using the level set method to evolve the interface Once we have obtained a discretization of the interface problem in space with second order accuracy, we have to move the interface. Among the several possibilities available, we have chosen the level set method [21]. Our problem is composed only by two phases, but similar procedures can be applied to multiphase flows [22]. Suppose that our problem presents a closed interface in the solution domain . An example can be the one shown in Fig. 1. For the simulation of the Hele–Show flow, this interface separates a flow with high viscosity from other flow with low viscosity. In the problems we present, the more viscous fluid will be inside the interface, i.e. in the “−” side. As shown in Fig. 1, we will call  (t) to the moving interface. − (t) and + (t) will be respectively the interior and exteriors regions of the interface.  (t) can be described as the zero isocontour of the level set function ϕ (r, z, t). This function is Lipschitz continuous and satisfies

ϕ(r, z, t) > 0

if (r, z) ∈ − ,

(69)

ϕ(r, z, t) = 0

if (r, z) ∈  ,

(70)

ϕ(r, z, t) < 0

if (r, z) ∈ + .

(71)

Differentiating the level set function in terms of t and taking into account that r and z will also depend on t, we obtain the Hamilton–Jacobi equation for the motion of the level set and of the interface

ϕt + u · ∇ϕ = 0

(72)

Initially we will choose ϕ (r, z, 0) as a signed distance function from the interface, i.e. |ϕ | = 1. Solving Eq. (72) we are able to update the position of the interface ϕ (r, z, t)=0.

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

191

Fig. 6. Least squares interpolation for the velocity at an irregular grid point.

4.3. Velocity field near the interface In order to solve the Hamilton–Jacobi equation to update the position of the interface, we will need the velocity field in order to obtain the normal velocity

un = −β pn = −β∇ p · n, where n is the normal vector at the interface. At regular grid points, the gradient can be obtained from the standard central difference scheme. The pressure is piecewise smooth, but the normal velocity might present a jump discontinuity determined by the natural jump condition v presented in (2). Thus, at irregular grid points we can use the IIM jointly with the least squares interpolation to obtain the velocity near the interface. In this case we would use the interface relations to express the points of the stencil from the wrong side of the interface in terms of the right side. Then we would use the Taylor expansion around the irregular grid point and the least squares interpolation in order to obtain the derivatives of the pressure, just as shown in [17,18]. In this case, the right hand side of the system that we need to solve in order to obtain γ i in Eq. (16), is different. For the derivative in terms of r we will use





γ1 , . . . , γ6 = 0  L2 γ1 , . . . , γ6 = cos θ   L3 γ1 , . . . , γ6 = − sin θ   L4 γ1 , . . . , γ6 = 0   L5 γ1 , . . . , γ6 = 0   L6 γ1 , . . . , γ6 = 0, L1



and for the derivative in terms of z we will use





γ1 , . . . , γ6 = 0   L2 γ1 , . . . , γ6 = sin θ   L3 γ1 , . . . , γ6 = cos θ   L4 γ1 , . . . , γ6 = 0   L5 γ1 , . . . , γ6 = 0   L6 γ1 , . . . , γ6 = 0, L1

as in this case we only want to find the gradient of the pressure at every irregular grid point. The solution of this system of equations will provide us with the weights γ i needed to obtain p. Once obtained γ i , we can obtain the correction term and then the gradient. Multiplying by the normal vector, we obtain the normal velocity at the projection. Now we only need to find some regular grid points in the neighborhood of the irregular grid point. At those points we have already obtained the velocity using central differencing. This information can be used to interpolate the normal velocity at the irregular grid point. If we choose two regular grid points plus the projection, as shown in Fig. 6, then we can say that the normal velocity at the irregular grid point (ri , zj ) is

un (ri , zj ) = α0 un (r0 , z0 ) + α1 un (r1 , z1 ) + α2 un (r2 , z2 ), where α 0 , α 1 , α 2 are the solutions of the system

α0 + α1 + α2 = 1,

(73)

α0 (ri − r0 ) + α1 (ri − r1 ) + α2 (ri − r2 ) = 0,

(74)

α0 (zi − z0 ) + α1 (zi − z1 ) + α2 (zi − z2 ) = 0.

(75)

192

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197 Table 4 A grid refinement analysis for the two components of the gradient (ur , uz ) with a different jump ratio in β − and β + . The numerical results are second order accurate. For the sake of brevity, we present only the results for the L norm. n

vr

vz

β = 1000, β = 1

β = 1, β = 1000

β + = 1000, β − = 1

β + = 1, β − = 1000

||EN ||

Order

||EN ||

Order

||EN ||

Order

||EN ||

Order

2.3822e−01 1.2175e−01 4.5039e−02 1.0280e−02 2.2699e−03 8.1657e−04

– 0.9465 1.4185 2.1194 2.1729 1.4729

9.7600e−01 5.5319e−01 2.1669e−02 7.8505e−03 2.0131e−03 5.6856e−04

– 0.8006 4.6214 1.4565 1.9578 1.8215

2.1596e−01 8.3636e−02 1.8423e−02 8.0693e−03 1.9525e−03 7.5113e−04

– 2.5822 4.5398 2.2831 4.1329 2.5994

8.5089e−01 4.8186e−01 2.8043e−02 7.6122e−03 2.0767e−03 6.1827e−04

– 1.3377 2.1580 1.1843 2.0414 1.3762

+

32 64 128 256 512 1024

−

+

−

−1

−1 y=1.7084x +4.2193

y=1.6519x +3.6378

−3

−3

−4

−4 log(|e∞|)

−2

log(|e∞|)

−2

−5

−5

−6

−6

−7

−7

−8 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

−8 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

Fig. 7. Linear regression plot for the components of the gradient vr (left) and vz (right) when β + = 1000 and β − = 1. In order to obtain this plot, we have used the data shown in Table 4. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

This is a second order interpolation scheme. More information can be found in [20] and [2]. This strategy would lead to an accurate enough approximation of the normal velocity. Using the solution and initial interface proposed in (62), we can check the accuracy of the solution obtained using the process described before. Tables 4 show a grid refinement analysis for the two components of the gradient. Once the solution of the Poisson problem  · (β u) = f has been obtained numerically, we can check the accuracy of the gradient u through a grid refinement analysis. The gradient of (62) is u = (ur , uz ), where

ur (r, z) = and

uz (r, z) =



2 2re(r ) · sin(z) if 3r2 if

 (r2 ) e · cos(z) if 3z2 if

ρ ≤ 0.4 + 0.04 cos(6θ ), ρ > 0.4 + 0.04 cos(6θ ), ρ ≤ 0.4 + 0.04 cos(6θ ), ρ > 0.4 + 0.04 cos(6θ ).

(76)

(77)

In Table 4 we check the numerical gradient with the previous expressions. We can observe nearly second order accuracy. Figs. 7 and 8 show the linear regression plot for the data shown in Table 4. 4.4. Reinitialization of the level set function We have chosen the level set function as a signed distance function from the interface. This means that |ϕ | = 1. The Hamilton-Jacobi equation will move the interface at the correct speed, but at the cost of transforming it into a function that is not a distance function any more. A reinitialization process is necessary in order to keep ϕ as a distance function. In our experiments we perform two to three reinitialization steps after every three to five time iterations. More information about the reinitialization process can be found in [23]

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197 1

193

1 y=2.1989x +6.8503

0

0

−1

−1

−2

−2

−3

−3

∞

log(|e |)

log(|e∞|)

y=2.2532x +7.1026

−4

−4

−5

−5

−6

−6

−7

−7

−8 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

−8 −7

−6.5

−6

−5.5

−5 log(h)

−4.5

−4

−3.5

−3

Fig. 8. Linear regression plot for the components of the gradient vr (left) and vz (right) when β + = 1 and β − = 1000. In order to obtain this plot, we have used the data shown in Table 4. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

5. Numerical experiments of the Hele–Shaw flow We have performed several experiments with different initial interfaces, viscosities and surface tensions. The results agree with those we can find in the literature. In particular, we show here some experiments inspired by those presented in [24,25] and [16]. Nevertheless, our aim is not to reproduce the experiments shown in previous references but to check the performance of the algorithm presented. In our case, the interface is tracked using a level set function. Thus, we cannot use the conformal mapping shown in [24] and [16] in order to model the initial interface. In our approximation, the initial interface is a circle plus a perturbation. The perturbation has been modeled as a Gaussian function with different values of amplitude and standard deviation. It is important to remark the difficulties that must be confronted when simulating the Hele–Shaw flow, due to the inner instability of the problem. In [24,25] and [16] the spatial resolutions used are much higher than the ones used in our experiments. The maximum resolution used in our experiments is N = 310 points while the maximum resolution used in the mentioned papers is N = 4096 points. However, the quality of our results is good. Also, it is important to bear in mind that the approach used in these articles, based on integral equations, allow us to obtain a higher resolution at the interface than using the level-set method. We want to simulate the movement of a blob of viscous fluid, i.e oil, surrounded by a less viscous fluid, i.e water, driven through a sink. In this case, the fluid interface collapses forming a finger during its motion towards the sink. For the experiment shown in Fig. 9 we present, at various times, the results obtained when the initial interface is a circle and the sink is placed at the point (r, z) = (0, 0.1). The equation of the initial interface can be written as the zero level set of the surface,



ϕ(r, z) = 1 − r2 + z2 . It can be observed that we present a reconstruction of the interface in two and three dimensions. In this case, the interface b2 is not perturbed. The surface tension is τ = 0.002. As mentioned before, the viscosity of the fluid has the expression μ = 12 β.

Inside the interface β − = 1 and outside the interface we set β + = 100. The domain has been set to r  (0, 2] and z  [ − 2, 2] − {0}. As R1 = 2, the spatial resolution is set to h = 2R1 /N, where N = 280 points. This means that we have 280 points in the z direction and 140 in the r direction. In order to avoid the singularity at (r, z) = (0, 0), the grid is staggered, as it was mentioned previously. The result of this experiment reproduces the classical one introduced in [24, 25], and [16]. It can be seen that a finger is developed as time advances and how the interface moves faster as the finger approaches the sink. The simulation has been stopped when the finger actually reaches the sink. It can be seen that it touches the sink far before the fluid is completely sucked. Fig. 10 shows another experiment. In this case, the initial interface is modeled through a circle perturbed with a Gaussian function. The sink is placed at the origin. It is inspired by the one circle experiment of [24]. The surface tension and the viscosity for this and the following experiments is always set to τ = 0.02, β + = 100 and β − = 1. The domain has been set to r  (0, 2] and z  [−2, 2] − {0}. In this case, N = 300. The equation of the initial interface can be written as the zero level set of the

194

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

Fig. 9. Simulation of a Hele–Shaw bubble with a sink placed at r = 0.1. The surface tension in this case is τ = 2e − 3. The initial interface corresponds to a circle of radius R = 1 centered at the origin. In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

Fig. 10. Simulation of a Hele–Shaw bubble with a sink placed at the origin. The surface tension in this case is τ = 2e − 3. The perturbation corresponds to a Gaussian function of amplitude A = 0.4 and standard deviations γ = π3 . In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

surface,

√ ⎧ ϕ(r, z) = 1 − r2 + z2 +  ⎪ ⎪ ⎪ ⎪ ⎪ ⎪   ⎨ φ − π /2 2 = A exp −  ⎪ γ ⎪ ⎪ ⎪

z ⎪ ⎪ ⎩φ = arctan r

(78)

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

195

Fig. 11. Simulation of a Hele–Shaw bubble with a sink placed at the origin. The surface tension in this case is τ = 2e − 3. The perturbation corresponds to a Gaussian function of amplitude A = 0.4 and standard deviations γ = π7 . In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

Fig. 12. Simulation of a Hele–Shaw bubble with a sink placed at the origin. The surface tension in this case is τ = 2e − 3. The perturbation corresponds to a Gaussian function of amplitude A = 0.4 and standard deviations γ = π4 . In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

196

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

Fig. 13. Simulation of a Hele–Shaw bubble with a sink placed at the origin. The surface tension in this case is τ = 2e − 3. The initial interface is (1.05z)2 + r2 /(1.05)2 = 1. In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

Fig. 14. Simulation of a Hele–Shaw bubble with a sink placed at the origin. The surface tension in this case is τ = 1e − 5. The perturbation corresponds to a sinusoidal function of amplitude A = 0.001. In the graph to the left, the horizontal axis is r and the vertical is z. The graph to the right corresponds to the surface of revolution that corresponds to the inner level set presented in the graph to the left. This surface is plotted in the Cartesian space XYZ, with z being the vertical axis. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)

We can see that the perturbation is modeled as a Gaussian function with a maximum at φ = π2 , where φ = arctan( zr ). π. The amplitude of this Gaussian has been set to A = 0.4 and the standard deviation has been set to γ = 16 Figs. 11 and 12 show again an initial interface perturbed as shown in (78), but using a standard deviations of the Gaussian functions as γ = π7 and γ = π4 respectively. The amplitude of the Gaussian is set to A = 0.5 and A = 0.4, respectively. These two experiments are inspired by the experiment two identical circles shown in [24]. In both of them we can see that two fingers are developed. The shape of the final interface before the fingers reach the sink has a high dependence on the perturbation introduced on the original interface. Fig. 13 shows the evolution of the initial interface (1.05r)2 + z2 /(1.05)2 = 1. We can write this interface as the zero level set of the surface,

"

ϕ(r, z) = 1 − (1.05z)2 +

r 2 . 1.05

(79)

J. Ruiz Álvarez, Z. Li / Applied Mathematics and Computation 264 (2015) 179–197

197

In Fig. 13 we observe that two fingers are developed and eventually reach the sink. This experiment is also inspired by the last experiment presented in [24]. Finally, Fig. 14 shows an experiment where the initial interface is a circle plus a sinusoidal perturbation of very small amplitude. More specifically, the initial interface can be written as,

√ ⎧ ϕ(r, z) = 1 − r2 + z2 +  ⎪ ⎪ ⎪ ⎪ ⎪ ⎨  = A sin(φ) ⎪ ⎪ ⎪

 ⎪ ⎪ ⎩φ = arctan z r

(80)

The value of the amplitude of the perturbation is A = 0.001 and the surface tension is τ = 1e − 5. Even though the perturbation is very small, we can see its influence in the final result. 6. Conclusions In this paper, we have developed the immersed interface method in cylindrical coordinates for problems with axial symmetry. The accuracy of the algorithm has been tested. A Matlab package for solving the PDE presented in (1) using the IIM has also been developed for solving this kind of problems. Some applications to the Hele–Shaw flow have been presented. Our intentions are to expand the work presented here in order to be able to find solutions to Stokes and Navier–Stokes equations for cylindrical problems with axial symmetry. Acknowledgment The first author would like to thank the coauthor and his university for the hospitality received during the author’s visits. The second author was partially supported by the US AFSOR grant FA9550-09-1-0520, the US NSF grant DMS-0911434, and the NIH grant 5R01GM96195-2, and CNSF grants 10971102, 11371199, and 11471166, BK20141443. The authors would like to thank the referees for their comments and suggestions which have helped to improve the quality of this manuscript. References [1] R.J. Leveque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (4) (1994) 1019–1044 doi:10.1137/0731054. [2] Z. Li, K. Ito, The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains (Frontiers in Applied Mathematics), Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2006. [3] J.R. Álvarez, J. Chen, Z. Li, The IIM in polar coordinates and its application to electro capacitance tomography problems, Numer. Algor. 57 (3) (2011) 405–423. [4] Y. Li, S.A. Williams, A.T. Layton, A hybrid immersed interface method for driven Stokes flow in an elastic tube, Numer. Math. Theor. Methods Appl. 6 (4) (2013) 600–616. [5] H. Nganguia, Y.-N. Young, An immersed interface method for axisymmetric Stokes flow, Private communications, 2014. [6] A. Gray, Modern Differential Geometry of Curves and Surfaces with Mathematica, 1st edition, CRC Press, Inc., Boca Raton, FL, USA, 1996. [7] Z. Li, A fast iterative algorithm for elliptic interface problems, SIAM J. Numer. Anal 35 (1995) 230–254. [8] T.Y. Hou, J.S. Lowengrub, M.J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys (1994) 312–338. [9] P.G. Saffman, G. Taylor, The penetration of a fluid into a porous medium or Hele–Shaw cell containing a more viscous liquid, Proc. Roy. Soc. Lond. Ser. A: Math. Phys. Sci. 245 (1242) (1958) 312–329. [10] J.-D. Chen, Radial viscous fingering patterns in Hele–Shaw cells, Exp. Fluids 5 (6) (1987) 363–371. [11] D. Pihler-Puzovic, R. Perillat, M. Russell, A. Juel, M. Heil, Modelling the suppression of viscous fingering in elastic-walled Hele–Shaw cells, J. Fluid Mech. 731 (2013) 162–183. [12] R.E. Goldstein, A.I. Pesci, M.J. Shelley, Instabilities and singularities in Hele–Shaw flow, Phys. Fluids (1994–present) 10 (11) (1998) 2701–2723. [13] E. Lajeunesse, J. Martin, N. Rakotomalala, D. Salin, Y. Yortsos, Miscible displacement in a Hele–Shaw cell at high rates, J. Fluid Mech. 398 (1999) 299–319. [14] A. Vasil, I. Markina, On the geometry of Hele–Shaw flows with small surface tension, Interfaces Free Boundar. 5 (2003) 183–192. [15] H.D. Ceniceros, T.Y. HOU, The singular perturbation of surface tension in Hele–Shaw flows, J. Fluid Mech. 409 (2000) 251–272 doi:10.1017/S0022112099007703. [16] H.D. Ceniceros, H. Si, Computation of axisymmetric suction flow through porous media in the presence of surface tension, J. Comput. Phys. 165 (1) (2000) 237–260 doi:10.1006/jcph.2000.6613. [17] H.D. Ceniceros, T.Y. Hou, H. Si, Numerical study of Hele–Shaw flow with suction, Phys. Fluids (1994–present) 11 (9) (1999) 2471–2486 doi:http://dx.doi.org/10.1063/1.870112. [18] Q. Nie, F.-R. Tian, Singularities in Hele–Shaw flows driven by a multipole, SIAM J. Appl. Math. 62 (2) (2001) 385–406. [19] R. Almgren, Singularity formation in Hele–Shaw bubbles, Phys. Fluids (1994–present) 8 (2) (1996) 344–352. [20] T.Y. Hou, Z. Li, S. Osher, H. Zhao, A hybrid method for moving interface problems with application to the Hele–Shaw flow, J. Comput. Phys. 134 (2) (1997) 236–252 doi:10.1006/jcph.1997.5689. [21] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, J. Comput. Phys. 79 (1) (1988) 12–49 doi:10.1016/0021-9991(88)90002-2. [22] H.-K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphase motion, J. Comput. Phys. 127 (1) (1996) 179–195 doi:http://dx.doi.org/10.1006/jcph.1996.0167. [23] S. Osher, R.P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Applied mathematical science, Springer, New York, N.Y., 2003. [24] Q. Nie, F.R. Tian, Singularities in Hele–Shaw flows, SIAM J. Appl. Math. 58 (1) (1998) 34–54 doi:10.1137/S0036139996297924. [25] B. Gustafsson, A. Vasilev, Conformal and Potential Analysis in Hele–Shaw Cell, Advances in Mathematical Fluid Mechanics, Springer, Basel, 2006.