Hele-Shaw Flow Simulation with Interactive Control using Complex Barycentric Coordinates

(1)

Ladislav Kavan and Chris Wojtan (Editors)

Hele-Shaw Flow Simulation with Interactive Control using Complex Barycentric Coordinates

Aviv Segall, Orestis Vantzos and Mirela Ben-Chen Technion - Israel Institute of Technology

Figure 1: (left) A typical setup of the Hele-Shaw experiment with our simulation results. (right) One of the effects obtained by our simulation.

Abstract

Hele-Shaw flow describes the slow flow of a viscous liquid between two parallel plates separated by a small gap. In some configurations such a flow generates instabilities known asSaffman-Taylor fingers, which form intricate visual patterns. While these patterns have been an inspiration for artists, as well as thoroughly analyzed by mathematicians, efficiently simulating them remains challenging. The main difficulty involves efficiently computing aharmonicfunction on a time-varying planar domain, a problem which has been recently addressed in the shape deformation literature using a complex-variable formulation of generalized barycentric coordinates. We propose to leverage similar machinery, and show how the model equations for the Hele-Shaw flow can be formulated in this framework. This allows us to efficiently simulate the flow, while allowing interactive user control of the behavior of the fingers. We additionally show that complex barycentric coordinates are applicable to the exteriordomain, and use them to simulate two-phase flow, yielding a variety of interesting patterns.

1. Introduction

The interaction between fluids often leads to compelling visual phenomena, such as mixing and pattern formation. In this paper we are interested inviscous fingering, which are the patterns gener- ated at the unstable interface of a viscous liquid. Such patterns can arise when a liquid flows into a porous medium (e.g. sand), and are closely related to other pattern phenomena such as bacterial growth and snowflake formation. One option to experimentally study such fingering phenomena, is to inject air into a viscous liquid trapped between two parallel plates separated by a small gap (see Figure2), also known as aHele-Shaw cell[Saf86]. This setup allows experimental and mathematical analysis of the pattern formation, as the governing equations for the expanding air bubble are the same as those of other more complex flows yielding similar phenomena.

From the Computer Graphics perspective, such flows generate intricate patterns which have inspired artists [Hal13] and design-

ers [Ner12]. It would therefore be potentially useful to simulate such patterns numerically, and allow the user to control the finger formation, while preserving the physical behavior and appearance of the liquid. While a plethora of methods exist for numerically simulating this phenomenon in the Computational Fluid Dynam- ics literature, the vast majority requires copious amounts of computational resources, and are thus not amenable to user control at interactive rates. Furthermore, traditional fluid simulation methods from Computer Graphics, such as a full Navier-Stokes simulation, is unnecessarily computationally heavy: there is no need to simulate the full behavior of the fluid in the domain, since the fingering phenomena happen at the moving free boundary.

In the spirit of recent methods for fluid simulation using boundary tracking [KB14], we suggest aboundary integralformulation for this problem. Our main observation is that the problem formulation shares many properties with the problem ofplanar shape

c

2016 The Author(s)

Eurographics Proceedings c2016 The Eurographics Association.

(2)

deformation, where the behavior is prescribed by user constraints, rather than by the laws of physics. We therefore propose to leverage a reduced model successfully used for shape deformation, namely generalized barycentric coordinates, in order to parameterize the behavior of the flow. As Hele-Shaw flow is governed by a harmonic function, we usecomplex holomorphic barycentric coordinates, which simplify the derivation and analysis.

We show how to formulate the model equations using complex barycentric coordinates, which allows us to simulate the flow at interactive rates, and thus allows user control over the direction in which the fingers grow. By controlling the domain of injection, e.g.

by injecting from aline segmentinstead of a point, we further the artist’s control and enable the generation of a large variety of patterns. Finally, we show that complex holomorphic coordinates are applicable to theexteriorof a planar bounded domain, which allows us to simulate finger formation in the case of two liquids with different viscosities, as well as for multiply connected domains, which allows us to simulate obstacles.

1.1. Related Work

While fingering in Hele-Shaw cells has not been, to the best of our knowledge, simulated in Computer Graphics, the body of work dedicated to the experimental, analytical and numerical study of this phenomenon in the Computational Fluid Dynamics (CFD) literature is vast, and a complete review is beyond our scope. We therefore focus our literature overview on putting our work in con- text of existing schemes, by discussing the simulation of this phenomenon in other disciplines, simulation of related phenomena in graphics, and other applications in graphics which use similar tools.

Viscous fingering in Hele-Shaw cells. An excellent review on the problem of viscous fingering in two dimensions, including the Saffman-Taylor model equations, the formulation using complex analysis and conformal maps, as well as numerical experiments, appears in [BKL^∗86]. A more recent mathematical treatment of the problem using a complex analytic approach is given in [GV06]. Ex- perimental investigation of this problem continues to this date, including, e.g., analyzing the dependency of the emerging pattern on the viscosity ratio in two-phase flow [BRN15]. Numerical methods in the CFD literature are diverse, including boundary integral methods [LLL07], yet the main focus in such disciplines is long-time evolution and the emergence of limit shapes (see e.g. the largest simulation to date [LLFPM09]), as opposed to computation at interactive rates which is necessary for enabling user control. For a recent review of numerical methods for this problem in CFD, see the PhD thesis [Dal13] and references within. Finally, it is worth noting that while the Cauchy integral formula has been used before [Kha15] for this problem, the formulation there is quite different, as the integral there is computed numerically as opposed to our approach which usesanalyticintegrals on polygonal domains, leading to a more stable computation.

Simulation of related phenomena in Graphics. For a review of the simulation of the full Navier-Stokes equations in graphics we refer to [Bri15]. The simulation of viscous flow using reduced dimensional methods has been proposed for viscous threads [BAV^∗10], viscous sheets [BUAG12] and viscous thin films on curved surfaces [AVW^∗15], and gap coupled solids [QYF15].

See e.g. [TDF^∗15], and references within, for additional approaches to viscous fluid simulation. As opposed to these methods, we only need to simulate the behavior of the boundary curve of the fluid, and therefore face different challenges. Perhaps the phenomenon most related to our approach, is the simulation ofLapla- cian growthleading to fractal pattern formation, which is governed by similar equations. Such phenomena are exhibited for example by lichen growth, as were simulated in [Sum01,DGA04] using Diffusion Limited Aggregation. In [KSSL07], a dielectric break- down model was used for efficiently simulating lightning, whereas in [Kim06] a hybrid algorithm was used for simulating ice formation. While all these problems are related to ours, the formulation of Hele-Shaw flow requires the use of dedicated solutions, which are both efficient and user controllable.

Other applications using similar tools. Our numerical simulation is based on complex-valued holomorphic barycentric coordinates, knowns as theCauchy-Green (CG) barycentric coordinates, which were first suggested for image deformation in [WBCG09], following their initial introduction using a real-variable formulation [LLCO08]. These coordinates were later extended to allow for conformal maps with sharp bends [WG10], to non-holomorphic functions [WBCGH11], and to three-dimensions [BCWG09]. The CG coordinates are a special case ofgeneralized barycentric coordinates, which are used in graphics mostly for cage-based shape deformation, see e.g. [Flo15], for a recent review. The CG coordinates are based on aboundary integral formulation, formulated in complex variables for ease of analysis, using analytical, as opposed to numerical, integration. It has been shown [WBCG09] that these coordinates are well-behaved even near the boundary of the domain, as they have a non-singular limit there, which motivates their use for the simulation of Hele-Shaw flows. Recently, beyond shape deformation, boundary element formulations have been used in graphics for, e.g., fluid simulation [KB14,BKB12,GNS^∗12], sound simulation [ZJ09], and crack simulation [HW15,ZBG15].

1.2. Contributions

Our main contribution is a formulation for efficiently simulating Hele-Shaw flow with viscous fingering at interactive rates, while allowing for user control, using Cauchy-Green barycentric coordinates. Specifically, we:

• Formulate the model equations of the Hele-Shaw flow in terms of the Cauchy-Green coordinates, which leads to an efficient numerical simulation method (Sections2,3).

• Show that the Cauchy-Green coordinates are applicable to more general problems, such as exterior domains, and multiply connected domains, which allows us to simulate two-phase flow, and flow with obstacles (Section4).

• Show a variety of effects that can be achieved with our technique (Section5).

2. One phase Hele-Shaw Flow 2.1. The Model.

The Physics. We investigate the evolution of an incompressible viscous liquid slowly injected into (or pumped out of) two parallel plates separated by a small gap, under the influence of surface

(3)

tension, and without gravity. To simplify the exposition, we ini- tially assume that the surrounding fluid is air (i.e. has zero viscosity and constant pressure), and extend later to more general settings.

We further assume no-slip boundary conditions at the interface between the liquid and the plates, and a freely evolving liquid-air interface. Figure2(a) illustrates this scenario.

The general Navier-Stokes equations describing fluid motion can be considerably simplified under the aforementioned assumptions.

Specifically, the fluid velocity can be integrated across the gap, yielding a reduced model in terms of the two-dimensional averaged velocityV. Following the derivation presented in [GV06], the governing equation isV =−∇Φ, whereΦis a scalarpotential function, related to the physical pressurepbyΦ= (^h²/12µ)p, with hthe gap height andµthe fluid viscosity.

Assuming the fluid is incompressible and fills the entire gap (therefore having a constant heighth) the fluid averaged velocity is divergence free everywhere except at the injection point, which we assume to be at the origin. There we have a source of strengthQ<0 representing a constant rate of injection. If the fluid is pumped out of the cell,Qwill be positive instead. Thus, in the interior of the fluid domain we have:

∆Φ=Qδ0(x,y),

where∆is the Laplacian andδ0(x,y)is the two-dimensional Dirac distribution supported at the origin.

The boundary conditions for the pressure p are given by the Young-Laplace condition, namely the pressure difference at the fluid-air interface is proportional to the mean curvature of the interface. Assuming constant air pressure at the exterior of the fluid, we can eliminate it from the equation by shifting both pressures by a constant factor. Furthermore, in the reduced two-dimensional model, the mean curvature of the interface is the curvatureκof the boundary curve, yielding the boundary conditionsΦ=σκ, where σis a rescaled surface tension parameter.

The Geometry. From the geometric perspective, the fluid occupies a time-dependent planar domainΩ(t)⊂C, which we assume to be simply-connected. Note, that we switch to complex-variable notation for points in thexy plane, namely we denote the point (x,y)∈R²byz=x+iy, whereiis the imaginary unit. The aforementioned model equations for the potential and velocity can be formulated as an evolution problem for the boundary of the domain

h Q<0 Q>0

(a)

Ω

𝑂

𝑛

𝜕Ω = Γ

(b)

Figure 2: The Hele-Shaw cell model. (a) The physical model. (b) The geometry and notation.

(a) (b) (c)

Figure 3: (top) injection and (bottom) suction, with zero surface tension. (a) The potential Φ at t =0 is positive for injection and negative for suction. (b) Boundary velocity: points closer to the source have higher velocity. (c) Curve evolution: the curve is smoothed for injection and sharpened for suction. The original curve is shown in blue, and later iterations in green.

Γ(t) =∂Ω(t), given in terms of the time-varying scalar potential Φ(t):Ω(t)→R[GV06, pp. 17]:

∆Φ(z) =Qδ0(z), z∈Ω (1a) Φ(z) =σκ(z), z∈Γ (1b) vn=h∂Γ

∂t(z),n(z)iˆ =h−∇Φ(z),n(z)i,ˆ z∈Γ, (1c) where ˆnis the outward unit normal direction of the boundary curve Γ(see Figure2(b)). The first two equations yield a unique solution for the potentialΦ(t), and the last equation specifies that the fluid- air interface (namely the boundaryΓ(t)of the domain) evolves according to the normal velocityvn=hV,ni.ˆ

Given an input initial domainΩ(0), our goal is to efficiently find a family of domainsΩ(t)which fulfill the model equations (1a)- (1c). To understand the behavior of the flow, consider the equations for thezero surface tensioncase (ZST), whenσ=0. In this case, the value ofΦon the boundary is 0, thus when the fluid is injected (i.e.Q<0), the potential in all the domain ispositive. Hence, the velocity at the boundary pointsoutwardand the boundary expands.

Intuitively, points closer to the singular point at the origin will have a larger potential gradient, and therefore move fasterawayfrom the origin. This effect tends tosmooththe curve. See Figure3(top) for an example showing the potential (a), the resulting velocity (b), and a few evolutions of the front under injection (c).

If, on the other hand, the fluid is pumped out (i.e.Q>0), the potential is negative in all the domain, and the velocity points to- wards theinterior. In this case as well points closer to the origin will move faster, but now the movement istowardsthe origin, en- hancing the curvature (see Figure3(bottom)). This property makes the front unstable, as small perturbations grow, and is the cause for the fingering phenomena. Numerically, this is one of the rea- sons simulating this flow is challenging: a naive discretization of the model equations in the case of suction (which is the interesting case generating the pleasing visual phenomena) might quickly be- come unstable and cease to evolve. While the surface tension term

c

2016 The Author(s)

(4)

acts as a regularizer, careful numerical treatment is still required in order to evolve the front in a stable and efficient manner.

To do that, we leverage an important property of the system, namely that it is described byharmonicfunctions, which allows us to reformulate the problem in terms ofboundaryinformation only. Specifically, we will consider two approaches, modeling the behavior ofΓandΦ, respectively. In both cases, reformulating the problem in terms of complex functions is instrumental, due to the wide applicability of complex methods to the analysis of harmonic problems in two-dimensions [CSMP15].

The Complex Formulation. We briefly mention some complex analysis notation which is required for the following discussion, and refer the reader to the excellent book [Ahl66] for a thorough introduction. We slightly abuse notation, by treating planar vectors (x,y)as the complex numberx+iy, thus for example, the gradient of a real functionφ:C→Rcorresponds to the complex number

∂φ/∂x+i^∂φ/∂y. Aholomorphicfunction is a function that iscomplex differentiable, namely the limit^∂^f/∂z(z₀) =limz→z0f(z)−f(z₀)/z−z₀

exists regardless to the direction in whichzapproachesz₀. TheCauchy-Riemannequations [Ahl66] formalize the relation between a holomorphic function f(z) =φ(z) +iψ(z)and its real and imaginary partsφ,ψ:C→R. Specifically,φ,ψare harmonic, and their gradients are orthogonal and of equal norm. Furthermore, any harmonic function is the real part of some holomorphic function. Thus, we can rephrase the Hele-Shaw model equations using a holomorphic complex potentialW:Ω→C, whose real part agrees with the real-valued potential: Re(W) =Φ. Reformulating Equa- tions (1a)-(1c) usingWwe have [GV06, pp.17-18]:

W(z) = Q

2πlog(z) +g(z), z∈Ω (2a)

Re(W(z)) =σκ(z), z∈Γ (2b)

vn=−Re(∂W

∂z n(z)),ˆ z∈Γ, (2c) wheregis a holomorphicregular function (i.e. without poles in Ω). For the first equation we used the fact that Re(¹/2πlog(z)) =

1/2πlog(|z|)is the Green’s function for the Laplacian in the plane, and thus solves Equation (1a), whereas g is used to fulfill the

𝑧 = 𝑓(𝑡, 𝜁)

𝑈(𝜁) Ω(𝑡, 𝑧)

𝑂

𝑛 =𝜁𝑓′ 𝜁 𝑛 = 𝜁 𝑓′ 𝜁

𝑂

𝜕Ω = Γ(t, z)

Figure 4: Notation for evolving the interface. We map the unit disk U(ζ)(left) using a time varying conformal map f(ζ,t)to a time- varying domainΩ(t,z)with boundaryΓ(t,z)(right). The normal to the disk is mapped with the derivative of the map f⁰to the scaled normal at the target domain.

analytic

(a)

ours

(b)

analytic

(c)

ours

(d) Figure 5: Comparison of the quadratic form analytic approach for injection (a) and suction (c) with our approach for injection (b) and suction (d), using the same initial curveΩ(0). Note that our method indeed produces a cusp similar to the cusp of the analytic solution.

boundary conditions (1b). Finally, the third equation is due to the representation of the inner product of two planar vectors in complex form:ha,bi=Re(ab), and the relation between the derivative of a holomorphic function and the gradient of its real part, yielding:

∂W/∂z=^∂Φ/∂x−i^∂Φ/∂y.

With the complex formulation at hand, we can now attempt to address the model equations. We will propose two approaches, with complementary advantages. First, we will leverage the invariance of harmonic functions underconformal (angle preserving) maps, to directly solve for the evolving frontΓ(t)by parameterizing it as a time-evolving conformal map (and thus a holomorphic function) from the unit disk. This allows us to handle both injection and suction, and produces similar behavior as a known analytic solution for the challenging case of suction with zero surface tension.

Unfortunately, this approach is difficult to extend to more general scenarios (e.g. non-zero surface tension and two-phase flow), and causes additional technical problems due to uneven sampling of the evolving front. Our second approach is to directly solve for the evolving complex potentialW, and it can be applied in a variety of scenarios, yet cannot reproduce the analytic solution. Still, this approach is highly useful in practice, as it is easily modified to allow for user control, and is efficient enough to allow interactivity.

Note that Figure5was produced with the first approach, and all the others were produced with the second approach.

2.2. Evolving the Interface

The Riemann Mapping Theorem states that for any simply connected domainΩ⊂Cthere exists a unique bijective conformal mapping which maps the unit diskU={ζ:|ζ|<1}intoΩsuch thatf:U→Ω,f(0) =0,f⁰(0)∈R⁺. Thus, we can track the time- varying domain of the fluidΩ(t)by the time-varying conformal map f(ζ,t)from the unit disk intoΩ(t)for everyt(see Figure4).

ThePolubarinova-Galin (PG) equation[GV06] provides a condition that the conformal mapping f(ζ,t)must satisfy (in the case of a single singular point at the origin (s=0) and zero surface tension) for the model equations to hold. It builds on three facts: First, harmonic functions are invariant under conformal maps, and thus given a solution to the model equations onUwe can use the conformal map to get a solution onΩ. Second, the normal velocityvn

can be expressed both in terms of the complex potentialW_Ωand the time derivative of the conformal map^∂^f/∂t. And finally, the normal

ˆ

nonΩcan also be computed usingf(as seen in Figure4).

Combining these facts yields the equation (see supplemental ma-

(5)

terial) [GV06, Eq. (1.16)]:

Re ∂f

∂tζ∂f

∂ζ

!

=−Q

2π, ζ∈∂U. (3)

It was shown [Gus84] that in the case of injection under some assumptions on smoothness of∂Ω(0)there exists a unique solution f(ζ,t)satisfying the PG equation. It is also possible to find analytic solutions by using a special form forf(ζ,t)(i.e. expressing specific types of boundaries). For example, in [Gal45] the author chose the quadratic formf(ζ,t) =a₁(t)ζ+a₂(t)ζ²wherea₁(t)anda₂(t)are real coefficients. Substituting f(ζ,t)into (3) gives two equations which can be solved for the coefficientsa1,a₂at timet, yielding an explicit solution for the problem.

In the next section we discuss the spatial discretization using the Cauchy-Green barycentric coordinates for this formulation, and the resulting discrete equations. Figure5shows such solutions to the PG equation for injection and suction, using the quadratic form approach and our approach, using the same initial curveΩ(0). Note, that our method produces similar behavior to the analytic solution.

2.3. Evolving the Potential

As solving for the conformal map f has several issues, we alter- natively suggest to find the complex potentialW(z)which satisfies Equations (2a)-(2c). We do so by solving for the holomorphic functiong(z):Ω→C, which satisfies the boundary conditions:

Re(g(z)) =−^Q/2πlog|z|+σκ(z). Interestingly, holomorphic functions and conformal maps are equivalent, thus we can use the same ansatz for the spatial discretization, namely the discrete Cauchy- Green coordinates. Furthermore, this approach is more easily gen- eralizable to handle multiple singularities of different types.

Mutliple Singularities. In the physical model, extending to multiple singularities implies that instead of having a single source or sink of the velocity at the origin, there are multiple sources and sinks at locationss_k∈Ω, with strengthsQ_k. Thus Equation (1a) changes to∆Φ=∑kQ_kδ(z−s_k). Since Green’s functions can be

(a) (b) (c)

Figure 6: Simulating a sink localized on a line segment. (a) The scalar potentialΦ. (b) The velocity of the interface, note the larger region of high velocities. (c) The resulting evolution of the front.

Figure 7: A few frames from an interactive simulation, where the user modifies the singularity’s location in real-time. The resulting singularity path is shown on the left. Note how the path of the fingers is modified to “aim” for the location of the closest singularity.

superimposed, the corresponding contribution to the complex potential is∑k1/2πQ_klog(z−s_k) =∑kWs(z,s_k).

Similar reasoning allows us to add line singularities, namely sources and sinks which are localized on line segments. Given a line segmentl:s(t) =z₁+t(z₂−z₁), its contribution to the complex potential isWl(z,l) =¹/2πQl

R1

0log(z−s(t))dt(see the supplemental material for the closed form solution of this integral). Fig.6 shows the scalar potential and the velocities for a source localized on a line segment. Note that, compared to a point source, there is a larger neighborhood of points on the evolving curve with large velocities, yielding a more noticeable effect during the evolution.

Combining the contributions from all the singularities yields to the following modification to Equation (2a):

W(z) =

∑

kWs(z,s_k) +

∑

kWl(z,l_k) +g(z), (4) where{s_k}and{l_k}are the sets of point sources and line segments, respectively.

In the next section we show how the Cauchy-Green barycentric coordinates can be used for this formulation. Figure7shows an example of a flow where the point location of the singularity (i.e. the sources) changes during the flow, which allows fine control on the behavior of the fingers. Since the computation is done at interactive rates, the user can move this location interactively, yielding an intuitive tool for generating finger-like effects (see the video).

3. Discretization

In the previous section we described how the model equations of Hele-Shaw flow can be reduced to finding a time-varying holomorphic function, representing either a conformal map from the unit disk to the fluid domain, or the regular part of the complex potential of the fluid domain, under some constraints. This setup is remarkably similar to the setup common inplanar shape deformation, where we seek a deformation of the input shape which isde- tail preserving, under some user constraints. In [WBCG09] it was proposed to use the machinery of conformal maps for this problem, yielding exactly the same mathematical formulation as we have, namely, finding a time varying conformal map under some constraints. We now leverage that machinery to get a deformation which is conformal, yet driven additionally by the physical model, rather than exclusively by a human user.

c

2016 The Author(s)

(6)

3.1. Cauchy-Green Coordinates.

TheCauchy integral formula[Bel92] is a central result in complex analysis, expressing the fact that the values of any holomorphic function inside a domainΩcan be calculated by the following integral on the boundary ofΩ:

f(z) = 1 2πi

I

∂Ω

f(w)

w−zdw, z∈Ω. (5)

The Cauchy-Green Coordinates [WBCG09] are a discretization of the Cauchy integral. The domainΩis discretized using a polygon on which we store the function as values at the vertices{fj}ⁿ_j=1. The function f(w)is approximated on each edge by a linear interpolation between these values. Then, the integration on the edges can be calculated analytically, yielding a complex coefficientCj(z) for each fj. These complex coefficients are called the Cauchy- Green barycentric coordinates. Finally, the integral is approximated using the sum:

f(z) =

∑

ⁿj=1Cj(z)fj.

Similarly, the derivative offcan be approximated using the derivative ofCj(z):

f⁰(z) =

∑

ⁿj=1C⁰_j(z)fj=

∑

ⁿj=1Dj(z)fj.

We provide the expression for the Cauchy-Green coordinates and their gradients in the supplemental material. In the following we show how the CG coordinates can be used for evolving the interface and the complex potential.

3.2. Evolving the Interface

Spatial Discretization. We search for a time varying conformal mapf:U→Ω, which satisfies Equation (3). We discretize the unit circle using a regularn-sided polygon ˆU, and represent the conformal map usingnfunctions fj(t),j∈1..n,t∈R, corresponding to the vertices of the polygon. Then, the map of ˆUis:

f(ζ,t) =

∑

ⁿj=1Cj(ζ)fj(t), ζ∈U.ˆ (6)

SinceCj(ζ)are independent of f, the time derivative is given by:_∂t^∂ f(ζ,t) =∑jCj(ζ)_∂t^∂fj. Thus, the semi-discrete PG equation corresponding to Equation (3) is, forζ∈∂U:ˆ

Re _n

j=1

∑

C_j(ζ)∂fj

∂t

! ζ

n m=1

∑

Dm(ζ)fm

!!

=−Q 2π. (7)

Time (seconds)

0 0.5 1 1.5 2

Area

3 3.5 4 4.5 5 5.5

Suction Injection

Figure 8: Front evolution for the stable case of injection (left) and suction (middle) with small surface tension (10⁻⁵), using the complex potential approach (Section 3.3). The method yields linear evolution of the area as expected (right).

We additionally sample the regular polygon at the pointsS= {ζ_l} ∈∂U, which leads to the space-discrete system of ODEs:ˆ

Re

(C∂

∂t

fˆ)_l(Df)ˆ _l

=−Q

2π,∀l∈1..|S|, (8) whereC,Dare complex matrices with entriesC_{l j} =Cj(ζ_l)and Dlm=ζlDm(ζ_l), respectively, and ˆfis a vector with entriesfj(t).

Time Discretization. We use an explicit Euler scheme to integrate equation (8). Specifically, given ˆf^kat iterationk, we find a discrete approximation of_∂t^∂ fˆ^k, denoted by(∆fˆ)^k, by minimizing the error of an over-constrained set of linear equations derived from Equa- tion (8) sampled at 4npoints. Finally, we set ˆf^k+1= fˆ^k+∆t(∆fˆ)^k for a constant delta time∆t=0.001. Figure5(top) shows a comparison of our evolution for the case of injection with the classic solution obtained using the quadratic complex form, where we achieve similar behaviour.

Regularization. While this approach works for the injection problem, the suction problem requires additional regularization because of its ill-posed nature. The regularization we propose is the mini- mization of the second spatial derivative of∆fˆ, in order to keep the conformal map smooth. Thus, we add a regularization termλC⁽²⁾ to the linear equations, whereC⁽²⁾ are the second spatial derivatives of the Cauchy-Green coordinates in matrix form (provided in the supplemental material). Figure5shows our result with this regularization (where we usedλ=0.001) compared with the classic analytic solution. Note that we manage to achieve the characteristic cusp despite our use of regularization.

3.3. Evolving the Potential

Spatial Discretization. We search for a holomorphic function g(z):Ω(t)→C, given by equations (4) and (2b). We first discretize the input domainΩ(t)usingnsamples, to get the closed polygon Ω(t), and then we use again the Cauchy-Green coordinates to rep-ˆ resentg(z):

g(z) =

∑

ⁿj=1Cj(z)g_j, z∈Ω(t).ˆ (9) The boundary conditions (2b) yield the constraints:

Re

n

∑

j=1

Cj(z)g_j

!

=−Re(Wsrc(z)) +σκ(z), z∈∂Ω,ˆ whereWsrc(z)is the combined potential of all the sources and sinks, as given in Equation (4).

We sample the boundary of the discrete domain at the points S={zl} ∈∂Ω, which again leads to an over-constrained systemˆ of linear equations, which can be solved for ˆg, the complex vector with entriesgj. The spatial derivative of the complex potential

∂/∂zWis computed using the known derivative of the potential at the singularities andg⁰(z_l) =∑jDj(z_l)gj. Finally, from Equation (2c), the normal velocity is given byvn=−Re(^∂/∂zWn), where ˆˆ nis the averaged normal at the vertices of ˆΩ.

Time Discretization. We use explicit Euler integration of Equa- tion (2c), and advance the sampled locations zj using z^k+1_j = z^k_j+ (∆t)vn(z_j)n(zˆ _j).Since in this setup we can directly prescribe non-zero surface tensionσ, no regularization is required. We do,

(7)

however, resample the curve∂Ω(t)ˆ during the evolution, taking into account the curvature. See section5.1for the details, as well as for the computation method of the dynamic time-step∆t.

While it is possible to use a more advanced time integrator, we have observed that this approach is efficient and stable. Specifically, for a constant rate of injectionQ, the area of the domain should grow linearly. Figure 8shows the result of injection (left) and suction (middle) from a single source using the complex potential approach and the corresponding graph denoting the change in the area (right). Note that we get a linear change in the area, as expected.

4. Extensions to the Model

The setup we presented, namely: simulating the one-phase Hele- Shaw flow by evolving the complex potential with the Cauchy- Green coordinates, can be easily extended to more complicated physical setups. We first present the generalization toexterior flow, namely the fluid occupies an unbounded domain in the plane which is the complement of a simply connected curve, by showing how the Cauchy-Green coordinates can be modified to handle holomorphic functions on unbounded domains. Then, by combining interior and exterior flows, we addresstwo-phase flowby solving for two potential functions. Finally, we show how to handleobstaclesusing multiply connected domains and different boundary conditions.

4.1. One Phase Hele-Shaw Flow with a Bubble

Ω 𝑛

𝜕Ω = Γ

We consider the inner fluid to be air with zero viscosity (forming a bubble inside the outer fluid) and suction or injection of the external fluid from infinity (see inset figure). The flow is driven by the potential of the external fluidΦwhich is related to the fluid pressure by a constant scal- ing factor. Since the singular point is at infinity the potential should be harmonic everywhere but behave at infinity as [DM13]:

Φ(z)≈ −Q

2πlog|z|, as |z| → ∞.

We therefore represent the potential asΦ(z) =−^Q/2πlog|z|+ h(z)whereh(z)is a harmonic function which tends to a constant at infinity, and its gradient tends to zero. The corresponding model equations are therefore:

W(z) =−Q

2πlog(z) +g(z), z∈/Ω (10a)

Re(W(z)) =−σκ(z), z∈Γ (10b)

vn=−Re(∂W

∂zn(z)),ˆ z∈Γ (10c)

whereg(z)is a holomorphic function which satisfies lim

|z|→∞g(z) = const, and ˆnstill points outward of the curve. As previously, the boundary conditions are given by the Young-Laplace condition re- lating the pressure difference to the curvature of the boundary. The viscosity of the inner fluid is negligible compared to the viscosity of

the outer fluid, and thus the pressure of the inner fluid is assumed to be constant, leading to the boundary conditions in Eq. (10b). Sim- ilarly to the interior flow, we will represent the holomorphic functiong(z)using the Cauchy-Green coordinates, by slightly modify- ing them to handle exterior domains.

Exterior Cauchy-Green Coordinates. Given a bounded simply connected domainΩand a function f(z)which is holomorphic in the exterior ofΩsuch that limz→∞f(z) =cfor some constantc, the following holds [Kas05, pp. 140]:

1 2πi

Z

∂Ω

f(w) w−zdw=

(

c z∈Ω

c−f(z) z∈/Ω.

This result is sometimes known asCauchy’s integral formula for an unbounded domain. Thus, we can pick an arbitrary pointa∈Ω, and then the value of f(z)for a pointz∈/Ωis given by:

f(z) = 1 2πi

Z

∂Ω

f(w) w−adw− 1

2πi Z

∂Ω

f(w) w−zdw, and is independent of the choice ofa.

We discretize the domain using a polygon ˆΩ, and use the Cauchy-Green coordinates for discretizing the Cauchy integral:

f(z) =

n

∑

j=1

Cê_j(z)fj, Cê_j(z):=Cj(a)−Cj(z), z∈/Ω,ˆ wherea∈Ωˆ is arbitrary, andCê_j(z)is theexterior Cauchy-Green coordinatefor a vertex jof ˆΩ. This result indicates that the exterior coordinates can be expressed using the regular Cauchy-Green coordinates, and so do their derivative asDê_j(z) =−Dj(z).

Exterior Flow. Using the exterior coordinates, we can apply our previous ansatz and discretize Equations (10a)-(10c). Specifically, we assume that lim_|z|→∞g(z) =constand represent it byg(z) =

∑jC^e_j(z)g_j. As before, the discrete values g_j are calculated by solving the over-constrained linear system obtained by sampling the boundary and applying the boundary conditions Re(g(z)) =

−σκ(z) +^Q/2πlog|z|. Given the values ofgj, the velocity is calculated using the derivative of the exterior coordinates and the interface is advanced according to the normal velocity.

Figure9shows an example of using the exterior flow to “con- tinue” a real Hele-Shaw flow. Specifically, we extracted from a photograph by the artist Antony Hall [Hal13] the boundary curve of a real Hele-Shaw flow, and used it as the initial conditions of our simulation. The Figure shows the original photograph (left), and our “simulated” photograph (right), after allowing the front to evolve (the initial fluid has darker color). Note that the simulated front closely resembles the original photograph. Fig.1,12,15and the attached video show additional results using the exterior flow.

4.2. Two Phase Hele-Shaw Flow

In the general case, there are two fluids with non-zero viscositiesµ₁ andµ₂occupying the interior and exterior of the domain [How00].

Their flow is driven by two harmonic potentials which we denote byΦ1andΦ2for the inner and outer fluid, respectively. We again represent the potentials as the real parts of complex holomorphic potentialsW1 andW2. For simplicity we assume a single source

c

2016 The Author(s)

(8)

Figure 9: “Continuing” an experimental Hele-Shaw flow. (left) Photograph by Antony Hall. (right) Our evolution starting from the boundary curve of the photograph.

inside the inner fluid located at the origin. Since the fluids are incompressible, the total amount of material must not change and thus injection of material at some location must be compensated by removal of material from another. Therefore, the outer fluid will also have a singularity, and we assume it is at infinity.

The corresponding equations for this model are [How00]:

W1(z) =Q1

2πlog(z) +g(z), z∈Ω (11a)

W2(z) =−Q₂

2πlog(z) +h(z), z∈/Ω (11b)

µ1Re(W₁(z))−µ2Re(W₂(z)) =σκ z∈Γ (11c)

−vn=Re ∂W1

∂z n(z)ˆ

=Re ∂W2

∂z n(z)ˆ

z∈Γ (11d) whereQ1is the strength of the singularity at the origin,Q2is the strength of the singularity at infinity,g(z)is a holomorphic function insideΩand h(z)is a holomorphic function in the exterior ofΩ which tends to a constant at infinity. Note that in order to preserve the incompressibility of the fluids we must have thatQ₁=−Q₂ (the rate of injection matches the rate of pumping). The holomorphic functionsg(z)andh(z)are determined by the Young-Laplace boundary condition (11c), expressing the pressure jump across the interface, and the kinematic boundary condition (11d), stating that the normal velocities of the two fluids at the interface must be equal (as the fluids do not mix).

As before, we represent the holomorphic functionsg(z)andh(z) using the interior and exterior Cauchy-Green coordinatesg(z) =

∑jC_j(z)g_jandh(z) =∑jC^e_j(z)h_j, and discretize equations (11a)- (11d) in the same way. The coefficientsgjandhjare found as the solution of an overconstrained linear system obtained by sampling the boundary, and the values ofgjare used for calculating the normal velocity and advance the boundary of the curve.

Figure10shows examples of two-phase flows, for the stable case of injection whenµ1>µ2(bottom), and the unstable case of injection for two viscosity ratios^µ1/µ₂(top,middle). Note that the smaller viscosity ratio generates more intricate and thinner fingers, as expected [BRN15]. Furthermore, it is worth noting that in the extreme limits of the viscosity ratio the two previous cases are recovered.

Specifically, when^µ2/µ1→0 the flow behaves as the one phase flow of the inner fluid and when^µ¹/µ2→0 it behaves as the one phase flow with a bubble. Figure13and the attached video show additional results of the two phase flow.

4.3. Obstacles.

Obstacles are formulated using the no-penetration Neumann boundary conditions, i.e. the normal velocity of the interface along the obstacle should be zero. Here the fluid domain may be multiply- connected, and its boundary∂Ωis composed of a free boundary denoted byΓ1, and a part which is allowed to move only in the tan- gent direction (where the boundary is part of an obstacle), denoted byΓ2. Thus, the formulation is similar to the formulation of the regular Hele-Shaw flow, with the exception that now the boundary condition for the potential function onΓ2is Re(^∂W_∂zn) =ˆ 0.

Since obstacles form holes in the domain, the domain is now multiply-connected. Interestingly, the Cauchy integral formula holds in this case as well [Bel92], with the modification that the orientation of the interior boundaries should be opposite to those of the exterior boundaries. Thus, we can use the same discretization as before using the Cauchy-Green coordinates to represent the regular part of the complex potential, and add the boundary conditions:

Re

Q 2πz+

n

∑

j=1

Dj(z)gj

! ˆ n

!

=0 z∈Γ2. Note, that in this case D(z) discretizes the multiply connected Cauchy integral. Given the fluid interface with the Cauchy-Green coordinates D^Γ(z) and m holes with the Cauchy-Green coordinates{D^k(z)}^m_k=1the multiply connected coordinates are given by D(z) =D^Γ(z)−∑^mk=1D^k(z). Fig.11 shows suction from a line source in an interior flow with obstacles and Fig.15shows an external flow with obstacles.

5. Experimental Results 5.1. Implementation details.

User Interface. We implemented our method in MATLAB. The interface is represented as a polygon withnvertices, wherenmay change during the flow. The user draws a control polygon, which is then interpolated using a cubic spline and sampled atn=100 points for getting the initial polygonal interface. The user adds singularity

t

Figure 10: Two phase flow simulation. (top) Unstable injection with

µ1/µ2=µa=0.01. (middle) Unstable injection with^µ1/µ2=µb= 0.3. (bottom) Stable injection with^µ¹/µ2=2. Note that the lower viscosity ratioµa(top) generates thinner and more intricate fingers.

(9)

Figure 11: Flow with obstacles, suction from an interior segment.

points and line singularities and chooses their strengthQ. The user is free to move the singularity locations during the simulation, and thus change the direction which the fingers will follow. Using the line singularities the user can choose the path of a finger when it reaches the line (see Figure12).

Simulation. For each simulation frame the interface is sampled at 4npoints (each edge is sampled 4 times) on which the boundary conditions are applied to obtain 4nlinear equations. The calcula- tion of the coordinates for these samples can be easily parallelized, and is thus done on the GPU (on an NVIDIA GTX 980 card). To give a feel for the timings involved, calculating the coordinates of 4000 points takes 5 milliseconds. The system of linear equations is then solved by minimizing the least squares error, resulting in a vector of coefficients representing the potential. The normal velocity at each vertex is then calculated using the derivative of the coordinates, and the vertices are moved using an explicit Euler scheme with a dynamic time step, which is chosen according to the ratio of the edge length and the normal velocity∆t=min(^|eⁱ^|/vn). Finally, we fit a cubic spline which interpolates the new polygon and sample it according to the curvature (i.e. more samples in the highly curved regions). The number of sampled points is chosen dynami- cally according to a minimal edge limit and a limit on the number of points, where we used 0.02 and 1000, respectively.

Singular integrals at the boundary. The Cauchy-Green coordinates and their derivatives can be singular when evaluated at the boundary of the domain. The coordinates, though, have a non- singular limit, given in [WBCG09], which we use for our com- putations. The derivatives have a non-singular limit on the edges of the boundary polygon, yet are undefined at the vertices. Thus, we calculate^∂W/∂zat a point close to the vertex inside the domain. We chose to calculate the derivative at a point with distance of 10⁻³ from the vertex in the normal direction into the interior or the exterior of the domain, depending on where the complex potential is defined (the interior or the exterior flow). In the two phase case we can calculate the velocity from either the interior or the exterior

Figure 12: By prescribing a line singularity the user controls the path of the fingers, as they follow the line when it is reached.

Figure 13: Unstable injection from the origin. See text for details.

potentials. The normal of the vertices is calculated as a weighted average of the incident edges normals.

Degrees of freedom. Since the coordinates sum to one, their imaginary parts sum to zero. Thus, we have one degree of freedom which can be fixed by choosing the imaginary part of the first coordinate to be zero. In the two phase case we have three degrees of freedom: two of them are expressed as a constant addition to the imaginary parts of each of the potentials, and fixed similarly. The third is due to the Young-Laplace boundary condition, as it involves the difference between the two potentials. It is fixed by choosing the real part of the first coordinate of one of the potentials to be zero.

5.2. Limitations.

Our method has a few limitations. First, we do not handle topology changes, which sometimes may be required (e.g. merging fronts after passing an obstacle, or bubbles created due to self intersec- tions). In principle, topology changes can be handled using a more sophisticated tracking algorithm. Second, for exterior flow, if the front becomes very large, the computational cost becomes larger as we require many points to represent the front. We believe that a multi-resolution approach, e.g. using a multi-grid based method could alleviate this problem, but leave further investigation for future work.

5.3. Applications.

Visualizing the interior flow with a texture. In this experiment we used a texture to visualize the flow in the interior of the domain.

We simulated the two phase case, where the boundary of the mesh acts as the interface between the two fluids. After solving for the potentials, we used the potential of the interior domain for moving the interior vertices of the mesh as well as the boundary vertices.

After each iteration we resample the boundary and the interior of the mesh and interpolate the texture coordinates. In Fig.13we show the results for unstable injection.

Pumping from the medial axis. Here we have computed the medial axis of an input curve, and used it as a collection of line singularities from which we pump the fluid (see Fig.14(left)).

Figure 14: (Left) Pumping fluid from the medial axis of the boundary. (Right) Directing the fingers by moving a point singularity.

c

2016 The Author(s)

(10)

Figure 15: Exterior flow with multiple obstacles.

Controlling the fingers. Here we control the direction which the fingers follow by moving the suction point in exterior flow. In Fig.14(right) the user moves the suction point in the shown path, and the fingers follow this path as shown in the next images. Note that in the figure we show the caged air in gray and do not show the fluid (which occupies the exterior of the domain).

Obstacles. In this experiment we tested exterior flow with multiple obstacles. In Figure15the fingers are forced to pass through the obstacles as they move toward a line source placed at the bottom.

The full simulation is showed in the attached video. Note that in this simulation we show the air inside the domain in red and do not show the fluid occupying the exterior of the domain.

6. Conclusions and Future Work

We proposed a method for interactively simulating Hele-Shaw flows using complex holomorphic barycentric coordinates. We demonstrated a variety of flow scenarios, such as interior, exterior and two phase, and showed how to incorporate obstacles. In addition, we provided a few applications for generating interesting curve deformations, and appealing texture effects.

We believe that our suggested approach, leveraging methods from shape deformation to be used for fluid simulation with interactive control, is quite promising, and could be generalized to other scenarios. For example, there exists a generalization to three dimensions of the Cauchy-Green coordinates, which could be potentially useful for three dimensional Laplacian growth. In addition, it might be possible to extend the interface tracking method to equations which are conformal invariant other than Laplace’s equation [Baz04]. It could also be potentially possible to apply other complex barycentric coordinates for other types of two dimensional flows, and for other applications in fluid simulation.

Acknowledgments

The authors acknowledge ISF grant 699/12, Marie Curie CIG 303511, and the German-Israeli Foundation for Scientific Research and Development, Grant No: I-2378-407.6

References

[Ahl66] AHLFORSL. V.:Complex analysis. 1966.4

[AVW^∗15] AZENCOTO., VANTZOSO., WARDETZKYM., RUMPFM., BEN-CHENM.: Functional thin films on surfaces. InProc. of SCA (2015), ACM, pp. 137–146.2

[BAV^∗10] BERGOU M., AUDOLYB., VOUGAE., WARDETZKYM., GRINSPUNE.: Discrete viscous threads. 116.2

[Baz04] BAZANTM. Z.: Conformal mapping of some non-harmonic functions in transport theory. InProceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences(2004), vol. 460, The Royal Society, pp. 1433–1452.10

[BCWG09] BEN-CHENM., WEBERO., GOTSMANC.: Variational harmonic maps for space deformation. InACM Trans. on Graphics (TOG) (2009), vol. 28, ACM, p. 34.2

[Bel92] BELLS. R.:The Cauchy transform, potential theory and conformal mapping. CRC press, 1992.6,8

[BKB12] BROCHUT., KEELERT., BRIDSONR.: Linear-time smoke animation with vortex sheet meshes. InProc. of SCA(2012), Eurographics Association, pp. 87–95.2

[BKL^∗86] BENSIMON D., KADANOFFL. P., LIANG S., SHRAIMAN B. I., TANGC.: Viscous flows in two dimensions. 977.2

[Bri15] BRIDSONR.: Fluid simulation for computer graphics. CRC Press, 2015.2

[BRN15] BISCHOFBERGERI., RAMACHANDRANR., NAGELS. R.: An island of stability in a sea of fingers: emergent global features of the viscous-flow instability. 7428–7432.2,8

[BUAG12] BATTYC., URIBEA., AUDOLYB., GRINSPUNE.: Discrete viscous sheets. 113.2

[CSMP15] CUMMINGSL., SMITHS. L., MARTINP., PROTASB.: Mod- ern applications of complex variables: Modeling, theory and computation.4

[Dal13] DALLASTONM. C.: Mathematical models of bubble evolution in a hele-shaw cell.2

[DGA04] DESBENOITB., GALINE., AKKOUCHES.: Simulating and modeling lichen growth. InComputer Graphics Forum(2004), vol. 23, Wiley Online Library, pp. 341–350.2

[DM13] DALLASTONM. C., MCCUES.: An accurate numerical scheme for the contraction of a bubble in a hele–shaw cell. 309–326.7 [Flo15] FLOATERM. S.: Generalized barycentric coordinates and appli-

cations. 161–214.2

[Gal45] GALINL.: Unsteady filtration with a free surface. InDokl. Akad.

Nauk USSR(1945), vol. 47, pp. 246–249.5

[GNS^∗12] GOLAS A., NARAIN R., SEWALL J., KRAJCEVSKI P., DUBEYP., LINM.: Large-scale fluid simulation using velocity-vorticity domain decomposition. 148.2

[Gus84] GUSTAFSSONB.: On a differential equation arising in a hele shaw flow moving boundary problem. 251–268.5

[GV06] GUSTAFSSONB., VASILÂA ´˘ZEVA.: Conformal and potential analysis in Hele-Shaw cells. Springer Science & Business Media, 2006.

2,3,4,5

[Hal13] HALL A.: Hele-shaw cell problem 2013. http://www.

antonyhall.net/works01/heleshaw02.html, 2013.1,7 [How00] HOWISONS. D.: A note on the two-phase hele-shaw problem.

243–249.7,8

[HW15] HAHND., WOJTANC.: High-resolution brittle fracture simulation with boundary elements. 151.2

[Kas05] KASANAH.:Complex Variables: Theory and Applications. PHI Learning Pvt. Ltd., 2005.7

[KB14] KEELER T., BRIDSON R.: Ocean Waves Animation using Boundary Integral Equations and Explicit Mesh Tracking. InProc. of SCA(2014), Koltun V., Sifakis E., (Eds.).1,2

[Kha15] KHALIDA.:Free boundary problems in a Hele-Shaw cell. PhD thesis, UCL (University College London), 2015.2

[Kim06] KIMT. W.-H.: Physically-based simulation of ice formation.

PhD thesis, University of North Carolina at Chapel Hill, 2006.2 [KSSL07] KIMT., SEWALLJ., SUDA., LINM. C.: Fast simulation of

laplacian growth. 68–76.2

(11)

[LLCO08] LIPMANY., LEVIND., COHEN-ORD.: Green coordinates.

InACM Trans. on Graphics (TOG)(2008), vol. 27, ACM, p. 78.2 [LLFPM09] LI S., LOWENGRUB J. S., FONTANA J., PALFFY-

MUHORAYP.: Control of viscous fingering patterns in a radial hele-shaw cell. 174501.2

[LLL07] LIS., LOWENGRUBJ. S., LEOP. H.: A rescaling scheme with application to the long-time simulation of viscous fingering in a hele–

shaw cell. 554–567.2

[Ner12] NERVOUSSYSTEM: Laplacian growth 2d.

http://n-e-r-v-o-u-s.com/projects/albums/

laplacian-growth-2d/, 2012.1

[QYF15] QIUL., YUY., FEDKIWR.: On thin gaps between rigid bod- ies two-way coupled to incompressible flow. Journal of Computational Physics 292(2015), 1–29.2

[Saf86] SAFFMANP.: Viscous fingering in hele-shaw cells. 73–94.1 [Sum01] SUMNERR. W.: Pattern formation in lichen, 2001.2 [TDF^∗15] TAKAHASHIT., DOBASHI Y., FUJISHIROI., NISHITA T.,

LINM. C.: Implicit formulation for sph-based viscous fluids. InCom- puter Graphics Forum(2015), vol. 34, Wiley Online Library, pp. 493–

502.2

[WBCG09] WEBER O., BEN-CHEN M., GOTSMAN C.: Complex barycentric coordinates with applications to planar shape deformation.

InComputer Graphics Forum (2009), vol. 28, Wiley Online Library, pp. 587–597.2,5,6,9

[WBCGH11] WEBERO., BEN-CHENM., GOTSMANC., HORMANN K.: A complex view of barycentric mappings. InComputer Graphics Forum(2011), vol. 30, Wiley Online Library, pp. 1533–1542.2 [WG10] WEBERO., GOTSMANC.: Controllable conformal maps for

shape deformation and interpolation. InACM Trans. on Graphics (TOG) (2010), vol. 29, ACM, p. 78.2

[ZBG15] ZHUY., BRIDSONR., GREIFC.: Simulating rigid body fracture with surface meshes. 150.2

[ZJ09] ZHENGC., JAMESD. L.: Harmonic fluids. 37.2

c

2016 The Author(s)