- Email: [email protected]
Numerical simulations for the motion of soap bubbles using level set methods
Available online at www.sciencedirect.com
Computers & Fluids 37 (2008) 524–535 www.elsevier.com/locate/compfluid
Numerical simulations for the motion of soap bubbles using level set methods Myungjoo Kang b
a,*
, Barry Merriman b, Stanley Osher
c
a Department of Mathematics, Seoul National University, Seoul, Republic of Korea Department of Human Genetics, University of California, Los Angeles, CA 90095, United States c Department of Mathematics, University of California, Los Angeles, CA 90095, United States
Available online 25 July 2007
Abstract In this paper, we develop appropriate equations for bubble motions depending on curvature using the level set methodology. Our method consists of an appropriate finite difference scheme for solving our model equations and a level set approach for capturing the complicated motion between the bubbles. Results indicate that our models and the level set methods can handle complicated interfacial motions and topology changes, and that they can numerically simulate many of the physical features of bubble motions. 2007 Elsevier Ltd. All rights reserved.
1. Introduction In this paper, we will derive the proper model equations and numerical algorithms for soap bubbles. It is generally accepted that many researchers in the field of CFD have difficulties in dealing with singularities of fluid interfaces. To resolve these difficulties, we will use a level set approach to compute the motion of soap bubbles. We intend to capture the interface using a level set approach instead of explicitly tracking it. Front tracking methods usually require that we add or subtract points dynamically. In [2], an interesting front tracking method which does not explicitly reposition the points of interface was devised. The investigators cited good results, but it seems hard to implement in three space dimensions. Moreover, topological changes cause difficulties, as with all tracking methods. In [4], a level set formulation for moving interfaces with curvature dependent velocity was introduced. The level set function is typically a smooth function, denoted as /. The level set formulation eliminates the problem of reposition-
ing the points to a moving grid and is capable of capturing geometric properties of highly complicated boundaries including topological changes without explicitly tracking the interfaces. An application of the level set formulation was used in [1,6] for incompressible fluid flows. They found that it was best, at least close to the front, to keep / as the signed distance from the front to prevent the development of steep or flat gradients in /. This can be done by solving a simple initial value problem for / which leaves the front location unchanged. We will use the level set approach to solve the problem for motion of soap bubbles in 2D and 3D. This includes area preserving motion by velocity or by acceleration. 2. Equations of motions In our study, we shall consider area preserving motion for soap bubbles with curvature dependent velocity or acceleration. 2.1. Motion with curvature dependent velocity
*
Corresponding author. E-mail addresses: [email protected] (M. Kang), [email protected] (B. Merriman), [email protected] (S. Osher). 0045-7930/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.07.002
Consider the equation u ¼ rjn ½pn
ð1Þ
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
where u is the velocity, r is the surface tension, j is the mean curvature, n is the normal vector at the front and [p] is the jump of pressure across the interface. We need to find [p] to preserve the area. The change of area is given by Z Z u n ds ¼ ðrj ½pÞ ds ð2Þ C
525
is that it must conserve total momentum of the bubble, i.e. summing FMT over all surface elements should give zero. This assumes that all mass transfer is within the bubble and not the surrounding environment. Given these cautions, we proceed to derive the equations. Newton’s law dP/dt = F for an element of the surface is
C
where C is the front of a bubble. From this, [p] can be rewritten as R r C j ds ½p ¼ R ð3Þ ¼ rj ds C
dðmomentum of elementÞ dt
where j is the average curvature at the front. Now, our governing equation becomes
If the element has area jdAj, mass per unit l, and velocity u (note, u need not be normal to the surface), then
u ¼ rðj jÞn
momentum of element ¼ ðljdAjÞu
ð8Þ
force per area ¼ ½pn rkn þ f
ð9Þ
ð4Þ
This equation is valid both 2D and 3D. 2.2. Motion with curvature dependent acceleration In this section, we will derive the proper equations of motions for a soap bubble in 3-dimensions (which gives the 2-dimensional results as a special case). We derive the basic equation of bubble motion by applying Newton’s law dP/dt = F where P = particle momentum, F = force on particle to a small element of the bubble surface. Care is required when doing this, because F must include all sources of momentum to the element – in particular, the source due to mass transfer. Let ds be a small piece of the bubble surface, with mass dm, velocity u, momentum P = dmu. Then Newton’s law is dðdmuÞ ¼F dt
ð7Þ
where [p] is the pressure jump across the surface, which is chosen to preserve the volume, r is surface tension, k curvature, n the outward unit normal and f is any additional source of momentum per unit area, which includes the possible mass transport effects described above. We do not assume that f is normal to the surface. Thus Newton’s law becomes dðljdAjuÞ ¼ ð½pn rkn þ f ÞjdAj dt Isolating l
ð6Þ
where FMT is the momentum source due to the mass transport, and F contains the standard forces. In particular, if the element mass dm is not constant in time, then we must consider this additional term, since there must be mass transfer. However, even if dm were constant, we could hypothesize a mass transfer process that resulted in momentum transfer as well. Thus FMT depends on the details of the mass transfer processes that actually occur. FMT could also be used to include the effects of evaporation or condensation on bubble momentum. The form of FMT is an additional physical axiom, independent of the rest of the problem. The only restriction on the form of the term
du dt
ð10Þ
and dividing out jdAj yields
du 1 dðljdAjÞ ¼ ½pn rkn þ f u dt jdAj dt ¼ ½pn rkn þ f þ F
ð5Þ
and we must specify F. If ds does not exchange any mass with other parts of the bubble, then the only momentum sources are the standard forces of pressure, surface tension, gravity, etc. However, if ds does exchange mass with the rest of the bubble surface, then the momentum brought in or carried away from ds by the mass transport process must also be included in the momentum balance: dðdmuÞ ¼ F þ FMT dt
¼ force on the element ¼ ðforce per areaÞðarea of elementÞ
ð11Þ
where F simply stands for the terms indicated. F can be broken into two parts by expanding out the derivative 1 dðj dhAjÞ 1 dl F ¼ lu þ ð12Þ jdAj dt l dt In order to close this system, we need an expression for ðaÞ
1 dðjdAjÞ ; jdAj dt
ðbÞ ½p;
ðcÞl and
dl ; dt
ðdÞ f
ð13Þ
The first depends only the velocity field, the second will be derived by the volume conservation constraint, the third requires additional specification. In order to represent a physical bubble, the total mass (surface integral of l) should be constant, which is an added constraint on the possible form of (c). 2.2.1. Calculation of djdAj dt Let dA denotes the vector area of a surface element, with scalar area jdAj. We can derive several equivalent formulas for the fractional rate of change of jdAj.
526
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
1 djdAj ¼ u nj jdAj dt 1 djdAj 1 ¼ rs u jdAj dt jdAj djdAj ¼ r u n Du n dt
n
ð15Þ
From this, we can immediately derive
ð16Þ
1 djdAj 1 d 1 d ¼ n dA dA dA ¼ 2 jdAj dt dt jdAj dt jdAj
where $s is the surface divergence operator, which can be written generally as $ Æ u n Æ Du Æ n when u is extended to be defined off the surface. Note that formula (14) is only valid when the velocity is normal to the surface. Also, in the special case of normal velocity u = unn, formula (16) reduces to (14). We will prove (15) and (16). Let dA be the rate of change of the vector area element. Let ds1 and ds2 by two-orthogonal tangent vectors whose cross product is dA, dA ¼ ds1 ds2 ¼ jds1 jjds2 jn
ð17Þ
Then we have dðdAÞ dðds1 Þ dðds2 Þ ¼ ds2 þ ds1 dt dt dt
ð19Þ
ð20Þ
Let v1 and v2 be unit vectors in the ds1,ds2 direction, ds1 = jds1jv1, ds2 = jds2jv2, v1 · v2 = n. This yields dðdAÞ ¼ jds1 jjds2 jððDu v1 Þ v2 þ v1 ðDu v2 ÞÞ dt
1 djdAj ¼ r u n Du n jdAj dt
where H(/) is the Heaviside function, 1 if / P 0 H ð/Þ ¼ 0 if / < 0
ð21Þ
dðxÞ ¼
d H ðxÞ dx
ð32Þ
The level set evolution equation is
Therefore,
ð22Þ
for any vector w, and carry out the cross product explicitly, we find
ð23Þ
So the normal component is
ð24Þ
where $s is the surface divergence operator. If u is a 3D vector, note that rs u ¼ ðv1 Du v1 þ v2 Du v2 þ n Du nÞ n Du n ¼ r u n Du n
ð30Þ
ð31Þ
in the sense of distribution. Then Z dV ¼ H ð/Þt dX dt X
w ¼ ðw v1 Þv1 þ ðw v2 Þv2 þ ðw nÞn
dðdAÞ ¼ jdAjðv1 Du v1 þ v2 Du v2 Þ dt ¼ jdAjrs u
ð28Þ
X
/t þ u r/ ¼ 0
n
ð27Þ
2.2.2. Calculation of [p] We derive [p] from the constraint that volume remains constant. We will give this derivation using the level set representation. The bubble volume can be written as Z Volume : V ¼ H ð/ÞdX ð29Þ
Now, if we write out vector Du Æ v1, Du Æ v2 in terms of its v1, v2, n component using
dðdAÞ ¼ jdAjðv1 Du v1 þ v2 Du v2 Þn dt þ ðstuff in v1 ; v2 planeÞ
1 djdAj jdAj dt
We will also use d(/) where the Dirac delta function
where Du is the full derivative of u. Thus we obtain dðdAÞ ¼ ðDu ds1 Þ ds2 þ ds1 ðDu ds2 Þ dt
ð26Þ
So that using our above result yields
ð18Þ
For any tangent vector ds, we have dðdsÞ ¼ change in u over segment ds ¼ Du ds dt
dðdAÞ ¼ jdAjðr u n Du nÞ dt
ð14Þ
ð25Þ
where $ is the usual 3D divergence. Thus we get an expression only involving the normal of the surface, as desired
ð33Þ
½H ð/Þt ¼ H 0ð/Þ/t ¼ H 0 ð/Þðu r/Þ ¼ H 0 ð/Þjr/ju n This yields Z dV ¼ H 0 ð/Þjr/ju n dX dt X
ð34Þ
ð35Þ
j ¼ 0, then we must To preserve volume, if we have dV dt t¼0 have dtd ðdV Þ ¼ 0 dt Z d dV d H 0 ð/Þðui @ i /ÞdX ¼ dt dt dt X Z Z ¼ d0 ð/Þ/t ui @ i / dX H 0 ð/Þðui Þt @ i / dX X X |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} I II Z dð/Þui @ i /t dX ð36Þ X |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} III
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
1 dð/Þ u ru þ ð½pn rjn þ f þ FÞ r/ dX l X Z Z 1 dð/Þjr/jdX þ ðu ruÞ ndð/Þjr/jdX ¼ ½p l X ZX 1 þ ð37Þ ðrj f n F nÞdð/Þjr/jdX l Z X ðIIIÞ ¼ dð/Þui @ i /t dX Z ZX @ i ðdð/Þui /t Þ dX d0 ð/Þ@ i /ui /t dX ¼ X X Z dð/Þ@ i ui /t dX Z ZX ¼ d0 ð/Þ@ i /ui /t dX ðr uÞðu nÞdð/Þjr/jdX ðIIÞ ¼
Z
X
X
ð38Þ Therefore, Z Z d dV 1 dð/Þjr/jdX ðu ruÞ ¼ ½p dt dt X l X Z ndð/Þjr/jdX þ ðr uÞðu nÞdð/Þjr/jdX X Z 1 ðrj f n F nÞdð/Þjr/jdX ð39Þ X l From our volume preserving constraint, we can get [p]. R ½p ¼
X
ððr uÞðu nÞ ðu ruÞ n l1 ðrj f n F nÞÞdð/Þjr/jdX R 1 dð/Þjr/jdX X l
ð40Þ
2.2.3. Specification of l, dl and f dt To completely specify the equations of motion, l; dl (condt vective derivative of l) and f (additional momentum sources, including mass transfer effects) are required. No matter how these are specified, the motion will conserve volume and momentum, since the equations were derived from those principles. However, a physical bubble R motion should also conserve the mass of bubble (M ¼ surface ljdAj), unless we are modeling a process like evaporation or condensation, which can alter the mass. There are four cases for l and f. Each of these cases also has simple restriction to 2D. Case 1. l=constant on the surface and in time (e.g. l = 1, say)This does not conserve bubble mass, since mass will vary in proportion to the area, which varies in time. Thus, it does not correspond to a physical bubble motion. Since it is artificial, we may as well also take f = 0, for simplicity. This choice results in the simplest set of equation R ððn Du nÞu ðu ruÞÞ ndð/Þjr/jdX R ½p ¼ rj þ X dð/Þjr/jdX X ð41Þ
where j is the average curvature over the surface area.
527
Thus, the final equation of motion becomes ut þ u ru ¼ rðj jÞ uðr n Du nÞ Z 1 þ n ðn Du nÞðu nÞdð/Þjr/jdX A X Z 1 ð42Þ n ðu ruÞ ndð/Þjr/jdX A X Case 2. l = constant over the surface, but varies in time to M conserve mass(l ¼ AðtÞ and M = 1, say) This is the simplest motion that conserves mass, and corresponds to a real bubble with a fluid surface. So the mass is mobile enough to stay uniformly distributed. In this case there is mass transfer to a surface element during the motion. Since dtd ðljdAjÞ is not zero, and we need to postulate a form of the additional source of momentum f. We imagine that when mass transfers into/out of the element being considered, it takes its momentum with it. This suggests that the force FMT on the element ds, as described above, is proportional to the local rate of change of element mass d ðljdAjÞ. Also, since we are considering here an dt idealized mass transfer that is so rapid it keeps the mass density uniform on the surface, it in some sense samples and averages things over the whole surface, so that the mass undergoing transfer carries with it the average velocity of the entire bubble. This intuition leads one to postulate that the additional momentum transfer to the surface element ds due to the mass transfer is F MT ¼ u
d ðljdAjÞ dt
ð43Þ
where u is the average vector velocity of the bubble Z 1 ujdAj ð44Þ u¼ A MT and, since f ¼ FjdAj ,
1 d ðljdAjÞ jdAj dt 1 dl ¼ lu r u n Du n þ l dt
f ¼u
ð45Þ
In this case, we can again make the same simplifications as in Case 1 to get [p], Z 1 ½p ¼ rj þ 2 ðu uÞ nðr u n Du n A X Z 1 dl 1 dð/Þjr/jdX þ 2 ððu nÞðr uÞ þ l dt A X ðu ruÞ nÞdð/Þjr/jdX ð46Þ and the equation of motion is
528
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
lðut þ u ruÞ ¼ rðj jÞn þ lðu uÞ r u n Z 1 dl n Du n þ 2 ðu uÞ l dt A X 1 dl n r u n Du n þ l dt Z n dð/Þjr/jdX 2 ððu nÞðr uÞ A X ðu ruÞ nÞdð/Þjr/jdX ð47Þ Case 3. ljdAj is constant during the motion (i.e. d ðljdAjÞ ¼ 0).This also conserves mass; it corredt sponds to an elastic membrane, where the mass is not free to migrate. In this case, we envision no mass transfer effect, thus f = 0. We get, from d ðljdAjÞ ¼ 0, that dt 1 dl 1 djdAj ¼ ¼ ðr u n Du nÞ l dt jdAj dt
ð48Þ
Therefore, R ððr uÞðu nÞ ðu ruÞ n l1 rjÞdð/Þjr/jdX X R1 ½p ¼ dð/Þjr/jdX l ð49Þ
Case 4. l flows with some mass preserving velocity on the surface. This velocity would generally require additional equations to specify it. The local momentum source due to mass transfer would be f = $s(lul), the surface divergence of momentum flux associated with l, velocity ul.
dA ! dL jdAj ! jdLj A!L Note that dL = njdLj is an outward normal. 3.3. Initial data In order to have volume preserving motions, the initial velocity of the bubble must also be volume preserving. ¼ 0attime ¼ 0). From Eq. (34), we see this means (i.e. dV dt that the initial velocity u0 must satisfy Z 0¼ u0 dA ð50Þ surface
If this condition is violated, there will be very large acceleration produced when the time evolution starts. The easiest way to insure (50) is to either take u0 = 0 which corresponds to a bubble at rest initially, or u0 Æ n = 0 which corresponds to a bubble ‘‘rotating’’ initially (i.e., velocity tangent to surface everywhere). 3.4. General equation of bubble motion The general equation for bubble motion is du 1 djdAj 1 dl þ l ¼ ½pn rjn þ f lu dt jdAj dt l dt
ð51Þ
where 1 jdAj d ¼ r u n Du n jdAj dt
ð52Þ
and R
3. Summary of bubble motion equations We summarize our equations of bubble motion for various cases. All equations are valid in both 2D and 3D. 3.1. Notational convention We review the notations used below: dA
the vector area of surface element in the direction normal to the surface jdAj scalar area of the surface element n the unit outward normal to the surface j mean curvature r surface tension dQ ¼ Q þ t R u rQ for any quantity Q dt R Qave ¼ A1 surface QjdAj where A is the total area (A ¼ jdAj) 3.2. 2D versions of equation The 2D versions of all equations are obtained by replacing area of bubble surface by length of the bubble curve in all relevant places:
½p ¼
X
ððr uÞðu nÞ ðu ruÞ n l1 ðrj f n F nÞÞdð/Þjr/jdX R 1 dð/Þjr/jdX Xl
ð53Þ
In addition, we specify conservation of mass Z Mass : M ¼ ljdAj ¼ constant in time
ð54Þ
and no self-induced force Z 0 ¼ f jdAj
ð55Þ
3.5. Volume preserving acceleration motion Consider du ¼ ½pn rjn dt
ð56Þ
This preserves the bubble volume if Z 1 ½p ¼ rj þ ððr uÞðu nÞ ðu ruÞ nÞdð/Þjr/jdX A X ð57Þ
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
This model is volume preserving, but it does not conserve momentum or mass for the bubble, so it is not a physical bubble motion. However it is useful to compare it with volume preserving velocity case. 3.6. Volume and momentum preserving model
529
This preserves bubble volume, mass and momentum, so it is a physically realistic bubble motion. It corresponds to a bubble modeled as an elastic membrane – i.e. no migration of mass within the surface. Various interesting effects can be created by starting with non-uniform mass distributions.
This is the model described in Section 2.2.3, Case 1, with l=1
4. Numerical method
du ¼ ½pn rjn uðr u n Du nÞ dt
4.1. Description of the level set approach
where ½p ¼ rj þ
R X
ð58Þ
ððn Du nÞðu nÞ ðu ruÞ nÞdð/Þjr/jdX R dð/Þjr/jdX X ð59Þ
This model does not conserve bubble mass, but it does conserve volume and momentum.
We construct a level set function / such that bubble interface is the zero level set of /. We also initialize / to be the signed distance from the interface such that / is positive inside the bubble and negative outside the bubble. This is easy to do using the method in [6]. We represent the interface by surface of bubble ¼ fxj/ðxÞ ¼ 0g /ðxÞ > 0 inside bubble
ð68Þ
/ðxÞ < 0 outside bubble
3.7. Fluid bubble This is the model described in Section 2.2.3, Case 2, with bubble mass M = 1 du 1 dl l ¼ ½pn rjn þ f lu r u n Du n þ ð60Þ dt l dt 1 ð61Þ lðtÞ ¼ AðtÞ 1 dl f ¼ lu r u n Du n þ ð62Þ l dt Z 1 ½p ¼ rj þ 2 ðu uÞ n r u n Du n A X Z 1 dl 1 þ dð/Þjr/jdX þ 2 ððu nÞðr uÞ l dt A X ðu ruÞnÞdð/Þjr/jdX ð63Þ This preserves bubble volume, mass and momentum and so is physically realistic bubble motion. It corresponds to a bubble with an ideal fluid surface on which mass is very mobile and maintains a uniform mass distribution. 3.8. Elastic membrane bubble This is the model described in Section 2.2.3, Case 3 du ¼ ½pn rjn ð64Þ dt lt¼0 ¼ a positive valued function on surface ð65Þ 1 dl ¼ ðr u n Du nÞ ð66Þ l dt R ððr uÞðu nÞ ðu ruÞ n l1 rjÞdð/Þjr/jdX X R1 ½p ¼ dð/Þjr/jdX l l
ð67Þ
The idea of the level set method is to move / with the correct speed u at the front using the following differential equation: /t þ u r/ ¼ 0
ð69Þ
Next, we reinitialize, using the algorithm [6]. keeping / to be signed distance, at least near the front. Additionally, we save computational time performing these calculations only near the front. There are several localization algorithms available; we use the relatively simple algorithm developed in [7]. 4.2. Notations on using level set function We can rewrite the variables by using the level set function /. r/ jr/j r/ j ¼ r jr/j RR jdð/Þjr/jdx dy j¼ RR dð/Þjr/jdx dy 8 > <1 n¼
ð70Þ ð71Þ ð72Þ if / > a
0 if / < a px > :1 / 1 otherwise 1 þ þ sin a p a ( 2 1 if j/j < a 1 þ cos p/ a 2a
H a ð/Þ ¼
da ð/Þ ¼ Z Zsurface volume
0
QjdAj ¼
Z
otherwise
Qdð/Þjr/jdX ZX Q dV ¼ QH a ð/ÞdV X
ð73Þ
ð74Þ ð75Þ ð76Þ
530
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
Osher and Sethian [4]. Consider the following differential equation:
where a is the prescribed ‘‘thickness’’ of the interface (usually 3 grid points in our calculations) and Q is any quantity.
/t þ u r/ ¼ 0
4.3. Outline of the numerical method
For u a given function of space and time, we obtain higher order accuracy by using ENO type schemes both in time and in space. For the time discretization, we use Heun’s method [3]. Let
We can now summarize our algorithm. Step 1. Initialize / (x, t) such that / is a signed distance function to the front. Step 2. Solve the governing equation and get the velocity u and update the level function /. The method for updating the level function / was proposed by 1
1
0
1
0 t = 0.0125
1
0 t = 0.025
1
1
0
0
0
1
1
0 t = 0.05
1
1
0
1
0 t = 0.0625
1
0 t = 0.1
1
1
1
0
0
0
0 t = 0.025
1
0 t = 0.05
1
1
1
0
0
0
0 t = 0.075
1
1
1
0
0
0
0 t=0
1
0 t = 0.025
1
1
1
1
0
0
0
1
0 t = 0.075
1
0 t = 0.1
1
1
1
1
0
0
0
1
0 t = 0.1
1
0 t = 0.125
0 t = 0.15
1
0 t = 0.175
1
1
1
1
0
0
0
1
0 t=0
1
0 t = 0.0125
1
1
1
1
0
0
0
1
0 t = 0.0375
1
0 t = 0.05
1
1
1
1
1
1
0
0
0
0
0
0
1
1
0 t = 0.125
1
0 t = 0.2
1
Fig. 3. Area preserving curvature dependent velocity: spiral relaxes towards circle.
1
0 t = 0.15
0 t = 0.05
0
0 t = 0.0875
1
1
1
Fig. 1. Area preserving curvature dependent velocity: circle is an equilibrium solution.
0 t=0
Then we apply Heun’s method.
1
0
0 t = 0.075
ð78Þ
0
1
0 t = 0.0375
o / ¼ L½/ ot
1
0
0 t=0
ð77Þ
0 t = 0.175
1
0 t = 0.2
1
Fig. 2. Area preserving curvature dependent velocity: ellipse relaxes to circle.
0 t = 0.075
1
0 t = 0.0875
1
0 t = 0.025
1
0 t = 0.0625
1
0 t = 0.1
1
Fig. 4. Area preserving curvature dependent velocity: square relaxes to circle.
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
/mþ1 ¼ /m DtL½/m 1 1 Dt /mþ1 ¼ /m þ /mþ1 L½/mþ1 ð79Þ 2 2 2 For the space discretization, first consider the motion of curvature dependent velocity.We can rewrite Eq. (79) as /t ðj jÞn r/ ¼ 0
When we apply the numerical scheme to a Hamilton–Jacobi equation, care is required. For the term corresponding jjr/j, we use ENO scheme(we will explain it at step 3). All the derivatives of the term jj$/j are approximated by central differences. The term jj$/j depends on the 1 multiplication and/or division by ð/2x þ /2y Þ2 , which might be close to zero near some points. This implies that we should be consistent with an approximation to j and j$/j in the term jj$/j to avoid unnecessary errors. For the curva-
ð80Þ
Using the relation (72), (80) becomes /t þ jjr/j ¼ jjr/j
531
ð81Þ
1
1
1
1
1
1
0
0
0
0
0
0
0 t=0
1
0 t = 0.00625
1
0 t = 0.0125
1
0 t=0
1
0 t = 0.00625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.01875
1
0 t = 0.025
1
0 t = 0.03125
1
0 t = 0.01875
1
0 t = 0.025
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.05
1
Fig. 5. Area preserving curvature dependent velocity: starfish relaxes towards circle.
0 t = 0.0375
1
0 t = 0.04375
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.00625
1
0 t = 0.0125
1
0 t=0
1
0 t = 0.00625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.01875
1
0 t = 0.025
1
0 t = 0.03125
1
0 t = 0.01875
1
0 t = 0.025
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.05
1
Fig. 6. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
1
0 t = 0.03125
1
0 t = 0.05
1
Fig. 7. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
1
0 t=0
0 t = 0.0125
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.0125
1
0 t = 0.03125
1
0 t = 0.05
1
Fig. 8. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
532
M. Kang et al. / Computers & Fluids 37 (2008) 524–535 Table 1 Ratio of area for final area to initial area
Fig. 1
Fig. 3
Fig. 5
Mesh size
32 · 32
64 · 64
128 · 128
Initial area Final area Final/initial Initial area Ffinal area Final/initial Initial area Final area Final/initial
0.7873 0.7788 0.9892 0.9762 0.9343 0.9570
0.7861 0.7815 0.9942 0.9747 0.9580 0.9829 0.9261 0.9241 0.9979
0.7856 0.7832 0.9969 0.9740 0.9688 0.9947 0.9263 0.9249 0.9985
Table 2 Ratio of area for final area to initial area
Fig. 6
Fig. 7
Fig. 8
Mesh size
128 · 128
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Ffinal area Final/initial
1.1000 1.0810 0.9827 0.7224 0.6955 0.9627 0.8806 0.8928 0.9424
Fig. 9. Three dimensional volume preserving curvature dependent velocity: symmetric pinch off, relaxes towards two spheres. 1
1
1
0
0
0
0 t=0
1
0 t = 0.0125
1
1
1
1
0
0
0
0 t = 0.0375
1
0 t = 0.05
1
1
1
1
0
0
0
0 t = 0.075
1
0 t = 0.0875
1
0 t = 0.025
1
0 t = 0.0625
1
0 t = 0.1
1
Fig. 11. Curvature dependent acceleration (volume preserving): circle remains circle.
Fig. 10. Three dimensional volume preserving curvature dependent velocity: non-symmetric pinch off, relaxes towards two spheres.
ture dependent acceleration, first we solve the governing equation to get the velocity u. After that we solve the Eq. (77) using the second-order ENO scheme[5] for the space. For the x-direction in 2D(or 3D),
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
8 ui;j ðDx /i;j þ h2 minmodðDx Dxþ /i;j ; Dx Dx /i;j Þ > > > < if u > 0 i;j u r/ ¼ > ui;j ðDxþ /i;j h2 minmodðDxþ Dxþ /i;j ; Dxþ Dx /i;j Þ > > : otherwise minmodðx; yÞ ¼
y-direction and the z-direction for 3D problems also. Step 3. Construct a new distance function /.Consider the following equation: d t ¼ signðd 0 Þð1 jrdjÞ with d 0 ¼ /n
ð82Þ
1
0
1
0
0 t=0
1
1
1
1
0
0
0
0
0 t = 0.125
1
0 t = 0.25
1
0 t=0
1
0 t = 0.0625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.375
1
0 t = 0.5
1
0 t = 0.625
1
0 t = 0.1875
1
0 t = 0.25
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.75
1
0 t = 0.875
1
0 t=1
1
Fig. 12. Curvature dependent acceleration (volume preserving): oscillating ellipse.
0 t = 0.375
1
0 t = 0.4375
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.125
1
0 t = 0.25
0 t=0
1
1
1
1
1
0
0
0
0
0 t = 0.375
1
0 t = 0.5
1
0 t = 0.625
1
0 t = 0.0125
1 1
0
0 t = 0.0375
1
1
1
0 t = 0.05
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.875
1
0 t=1
1
Fig. 13. Curvature dependent acceleration (volume preserving): oscillating square.
1
0 t = 0.3125
1
0 t = 0.5
1
0 t = 0.075
1
0 t = 0.025
1
0 t = 0.0625
1
0
1
0 t = 0.75
0 t = 0.125
Fig. 14. Curvature dependent acceleration (volume preserving): oscillating starfish.
1
0 t=0
ð84Þ
where sign(d0) gives the sign of d0. Given initial data, we solve the above equation until the solution reaches a steady state near the front. To eliminate the stiffness of sign function, we approximate sign(/) by
signðxÞ minðjxj; jyjÞ if x y > 0 ð83Þ 0 otherwise
and Dx ; Dxþ are just backward and forward divided differences. Similarly, we can apply this to the 1
533
0 t = 0.0875
1
0 t = 0.1
Fig. 15. Curvature dependent acceleration (volume and momentum preserving): circle remains circle.
534
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
/ signð/Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi / 2 þ 2
ð85Þ
where is very small number (i.e. = h). Following [3,6], we use the approximation: 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 > > maxððaþ Þ ; ðb Þ Þ þ maxððcþ Þ ; ðd Þ Þ > > > > > < if / > 0 0 jr/j ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > maxðða Þ2 ; ðbþ Þ2 Þ þ maxððc Þ2 ; ðd þ Þ2 Þ > > > > : otherwise ð86Þ h minmodðDx Dxþ /i;j ; Dx Dx /i;j Þ 2 h b ¼ Dxþ / minmodðDxþ Dxþ /i;j ; Dxþ Dx /i;j Þ 2
a ¼ Dx / þ
1
1
1
0
0
0
0 t=0
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.5
1
1
1
1
0
0
0
0 t = 0.75
1
0 t = 0.875
1
0 t = 0.25
1
0 t = 0.625
1
0 t=1
1
Fig. 17. Curvature dependent acceleration (volume and momentum preserving): oscillating square.
a ¼ minða; 0Þ bþ ¼ maxðb; 0Þ b ¼ minðb; 0Þ
1
1
1
0
0
0
ð88Þ
We can define c, d using same approximation in the y-direction. For the criterion for steady state, we use Q¼
0 t = 0.125
1
ð87Þ
aþ ¼ maxða; 0Þ
1
n Rjd ni;j j
M
< ðDtÞh2
0 t=0
1
0 t = 0.0625
1
1
1
1
0
0
0
0 t = 0.125
1
0 t = 0.3125
1
0 t = 0.5
1
ð89Þ jd ni;j j
where M = number of grid point where < að¼ chÞ, for some constant c (we use c = 1.5). This
1
1
1
0
0
0
0 t=0
1
0 t = 0.125
1
0 t = 0.25
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.5
1
1
1
0
0
0
0 t = 0.75
1
0 t = 0.875
1
1
0 t = 0.25
1
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.4375
1
Fig. 18. Curvature dependent acceleration (volume and momentum preserving): oscillating starfish.
is a very local criterion and convergence is very rapid.After solving this equation we let 0 t = 0.625
1
1
0 t = 0.1875
1
/nnew ¼ steady state solution of ð84Þ:
ð90Þ
Step 4. We have now advanced one time step. Go to step 2 and repeat.
0 t=1
1
Fig. 16. Curvature dependent acceleration (volume and momentum preserving): oscillating ellipse.
5. Numerical results in 2D and 3D Our experiments simulate the motions of soap bubbles. In the case of curvature dependent velocity and accelera-
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
tion, we simulate the merging and breaking of the interface in 2D and 3D. 5.1. Results using curvature dependent velocity We simulate several different shapes in the case of curvature dependent velocity. In Figs. 1–5, we show the continuous evolution of non-intersecting curves collapsing smoothly to a circle. In Figs. 6–8, we show the merging of two bubbles in 2D. After merging they converge to a steady state which is, of course, a circle. Table 3 Ratio of area for final area to initial area
Fig. 11
Fig. 12
Fig. 13
Fig. 14
Mesh size
32 · 32
64 · 64
128 · 128
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial
0.7873 0.7806 0.9915 0.8792 0.7939 0.9030 1.0000 0.9856 0.9856 0.8873 0.7964 0.8976
0.7861 0.7831 0.9962 0.8790 0.8413 0.9570 1.0000 1.0240 1.0240 0.9201 0.8470 0.9205
0.7856 0.7841 0.9980 0.8791 0.8635 0.9823 1.0000 1.0320 1.0320 0.9267 0.8882 0.9582
Table 4 Ratio of area for final area to initial area
Fig. 15
Fig. 16
Fig. 17
Fig. 18
Mesh size
32 · 32
64 · 64
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial
0.7873 0.7806 0.9915 0.8792 0.7927 0.9017 1.0000 0.9810 0.9810 0.8873 0.8105 0.9135
0.7861 0.7831 0.9962 0.8790 0.8414 0.9572 1.0000 1.0130 1.0130 0.9201 0.8744 0.9503
535
In Figs. 9 and 10, we use a symmetric and non-symmetric dumbbell as initial data. In 2D, a dumbbell under our motion collapses to a circle. However, in 3D, a dumbbell breaks into two part and each part changes smoothly to a sphere. Both 2D and 3D motions preserve the area and the volume, respectively (see Tables 1 and 2). 5.2. Results using curvature dependent acceleration In Figs. 11–14, we simulate the motion of bubble which was previously described in Section 2.2.3, Case 1. The result shows that motion has periodic oscillations as expected. We next use the motion described in Section 2.2.3, Case 2 (see Figs. 15–18). Both cases preserve area well (see Tables 3 and 4). Acknowledgements M. Kang was supported by Korea Research Foundation Grant (R08-2004-000-10233-0). M. Kang also holds joint appointment in the Research Institute of Mathematics, Seoul National University. S. Osher was supported by ONR Grant N00014-03-1-0071. References [1] Chang YC, Hou TY, Merriman B, Osher S. A level set formulation of Eulerian interface capturing methods for incompressible fluid flows. J Comput Phys 1996;124:449–64. [2] Hou TY, Lowengrub J, Shelley M. Removing the stiffness interfacial flows with surface tension. J Comput Phys 1994;114:312–38. [3] Osher S, Shu C-W. High-order essentially nonoscillatory schemes for Hamilton–Jacobi equation. SINUM 1991;28:907–22. [4] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed : algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49. [5] Shu C-W, Osher S. Efficient implementation of essentially nonoscillatory shock capturing schemes II. J Comput Phys 1989;83:32– 78. [6] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59. [7] Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE based fast local level set method. J Comput Phys 1999;55:410–38.
Computers & Fluids 37 (2008) 524–535 www.elsevier.com/locate/compfluid
Numerical simulations for the motion of soap bubbles using level set methods Myungjoo Kang b
a,*
, Barry Merriman b, Stanley Osher
c
a Department of Mathematics, Seoul National University, Seoul, Republic of Korea Department of Human Genetics, University of California, Los Angeles, CA 90095, United States c Department of Mathematics, University of California, Los Angeles, CA 90095, United States
Available online 25 July 2007
Abstract In this paper, we develop appropriate equations for bubble motions depending on curvature using the level set methodology. Our method consists of an appropriate finite difference scheme for solving our model equations and a level set approach for capturing the complicated motion between the bubbles. Results indicate that our models and the level set methods can handle complicated interfacial motions and topology changes, and that they can numerically simulate many of the physical features of bubble motions. 2007 Elsevier Ltd. All rights reserved.
1. Introduction In this paper, we will derive the proper model equations and numerical algorithms for soap bubbles. It is generally accepted that many researchers in the field of CFD have difficulties in dealing with singularities of fluid interfaces. To resolve these difficulties, we will use a level set approach to compute the motion of soap bubbles. We intend to capture the interface using a level set approach instead of explicitly tracking it. Front tracking methods usually require that we add or subtract points dynamically. In [2], an interesting front tracking method which does not explicitly reposition the points of interface was devised. The investigators cited good results, but it seems hard to implement in three space dimensions. Moreover, topological changes cause difficulties, as with all tracking methods. In [4], a level set formulation for moving interfaces with curvature dependent velocity was introduced. The level set function is typically a smooth function, denoted as /. The level set formulation eliminates the problem of reposition-
ing the points to a moving grid and is capable of capturing geometric properties of highly complicated boundaries including topological changes without explicitly tracking the interfaces. An application of the level set formulation was used in [1,6] for incompressible fluid flows. They found that it was best, at least close to the front, to keep / as the signed distance from the front to prevent the development of steep or flat gradients in /. This can be done by solving a simple initial value problem for / which leaves the front location unchanged. We will use the level set approach to solve the problem for motion of soap bubbles in 2D and 3D. This includes area preserving motion by velocity or by acceleration. 2. Equations of motions In our study, we shall consider area preserving motion for soap bubbles with curvature dependent velocity or acceleration. 2.1. Motion with curvature dependent velocity
*
Corresponding author. E-mail addresses: [email protected] (M. Kang), [email protected] (B. Merriman), [email protected] (S. Osher). 0045-7930/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.07.002
Consider the equation u ¼ rjn ½pn
ð1Þ
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
where u is the velocity, r is the surface tension, j is the mean curvature, n is the normal vector at the front and [p] is the jump of pressure across the interface. We need to find [p] to preserve the area. The change of area is given by Z Z u n ds ¼ ðrj ½pÞ ds ð2Þ C
525
is that it must conserve total momentum of the bubble, i.e. summing FMT over all surface elements should give zero. This assumes that all mass transfer is within the bubble and not the surrounding environment. Given these cautions, we proceed to derive the equations. Newton’s law dP/dt = F for an element of the surface is
C
where C is the front of a bubble. From this, [p] can be rewritten as R r C j ds ½p ¼ R ð3Þ ¼ rj ds C
dðmomentum of elementÞ dt
where j is the average curvature at the front. Now, our governing equation becomes
If the element has area jdAj, mass per unit l, and velocity u (note, u need not be normal to the surface), then
u ¼ rðj jÞn
momentum of element ¼ ðljdAjÞu
ð8Þ
force per area ¼ ½pn rkn þ f
ð9Þ
ð4Þ
This equation is valid both 2D and 3D. 2.2. Motion with curvature dependent acceleration In this section, we will derive the proper equations of motions for a soap bubble in 3-dimensions (which gives the 2-dimensional results as a special case). We derive the basic equation of bubble motion by applying Newton’s law dP/dt = F where P = particle momentum, F = force on particle to a small element of the bubble surface. Care is required when doing this, because F must include all sources of momentum to the element – in particular, the source due to mass transfer. Let ds be a small piece of the bubble surface, with mass dm, velocity u, momentum P = dmu. Then Newton’s law is dðdmuÞ ¼F dt
ð7Þ
where [p] is the pressure jump across the surface, which is chosen to preserve the volume, r is surface tension, k curvature, n the outward unit normal and f is any additional source of momentum per unit area, which includes the possible mass transport effects described above. We do not assume that f is normal to the surface. Thus Newton’s law becomes dðljdAjuÞ ¼ ð½pn rkn þ f ÞjdAj dt Isolating l
ð6Þ
where FMT is the momentum source due to the mass transport, and F contains the standard forces. In particular, if the element mass dm is not constant in time, then we must consider this additional term, since there must be mass transfer. However, even if dm were constant, we could hypothesize a mass transfer process that resulted in momentum transfer as well. Thus FMT depends on the details of the mass transfer processes that actually occur. FMT could also be used to include the effects of evaporation or condensation on bubble momentum. The form of FMT is an additional physical axiom, independent of the rest of the problem. The only restriction on the form of the term
du dt
ð10Þ
and dividing out jdAj yields
du 1 dðljdAjÞ ¼ ½pn rkn þ f u dt jdAj dt ¼ ½pn rkn þ f þ F
ð5Þ
and we must specify F. If ds does not exchange any mass with other parts of the bubble, then the only momentum sources are the standard forces of pressure, surface tension, gravity, etc. However, if ds does exchange mass with the rest of the bubble surface, then the momentum brought in or carried away from ds by the mass transport process must also be included in the momentum balance: dðdmuÞ ¼ F þ FMT dt
¼ force on the element ¼ ðforce per areaÞðarea of elementÞ
ð11Þ
where F simply stands for the terms indicated. F can be broken into two parts by expanding out the derivative 1 dðj dhAjÞ 1 dl F ¼ lu þ ð12Þ jdAj dt l dt In order to close this system, we need an expression for ðaÞ
1 dðjdAjÞ ; jdAj dt
ðbÞ ½p;
ðcÞl and
dl ; dt
ðdÞ f
ð13Þ
The first depends only the velocity field, the second will be derived by the volume conservation constraint, the third requires additional specification. In order to represent a physical bubble, the total mass (surface integral of l) should be constant, which is an added constraint on the possible form of (c). 2.2.1. Calculation of djdAj dt Let dA denotes the vector area of a surface element, with scalar area jdAj. We can derive several equivalent formulas for the fractional rate of change of jdAj.
526
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
1 djdAj ¼ u nj jdAj dt 1 djdAj 1 ¼ rs u jdAj dt jdAj djdAj ¼ r u n Du n dt
n
ð15Þ
From this, we can immediately derive
ð16Þ
1 djdAj 1 d 1 d ¼ n dA dA dA ¼ 2 jdAj dt dt jdAj dt jdAj
where $s is the surface divergence operator, which can be written generally as $ Æ u n Æ Du Æ n when u is extended to be defined off the surface. Note that formula (14) is only valid when the velocity is normal to the surface. Also, in the special case of normal velocity u = unn, formula (16) reduces to (14). We will prove (15) and (16). Let dA be the rate of change of the vector area element. Let ds1 and ds2 by two-orthogonal tangent vectors whose cross product is dA, dA ¼ ds1 ds2 ¼ jds1 jjds2 jn
ð17Þ
Then we have dðdAÞ dðds1 Þ dðds2 Þ ¼ ds2 þ ds1 dt dt dt
ð19Þ
ð20Þ
Let v1 and v2 be unit vectors in the ds1,ds2 direction, ds1 = jds1jv1, ds2 = jds2jv2, v1 · v2 = n. This yields dðdAÞ ¼ jds1 jjds2 jððDu v1 Þ v2 þ v1 ðDu v2 ÞÞ dt
1 djdAj ¼ r u n Du n jdAj dt
where H(/) is the Heaviside function, 1 if / P 0 H ð/Þ ¼ 0 if / < 0
ð21Þ
dðxÞ ¼
d H ðxÞ dx
ð32Þ
The level set evolution equation is
Therefore,
ð22Þ
for any vector w, and carry out the cross product explicitly, we find
ð23Þ
So the normal component is
ð24Þ
where $s is the surface divergence operator. If u is a 3D vector, note that rs u ¼ ðv1 Du v1 þ v2 Du v2 þ n Du nÞ n Du n ¼ r u n Du n
ð30Þ
ð31Þ
in the sense of distribution. Then Z dV ¼ H ð/Þt dX dt X
w ¼ ðw v1 Þv1 þ ðw v2 Þv2 þ ðw nÞn
dðdAÞ ¼ jdAjðv1 Du v1 þ v2 Du v2 Þ dt ¼ jdAjrs u
ð28Þ
X
/t þ u r/ ¼ 0
n
ð27Þ
2.2.2. Calculation of [p] We derive [p] from the constraint that volume remains constant. We will give this derivation using the level set representation. The bubble volume can be written as Z Volume : V ¼ H ð/ÞdX ð29Þ
Now, if we write out vector Du Æ v1, Du Æ v2 in terms of its v1, v2, n component using
dðdAÞ ¼ jdAjðv1 Du v1 þ v2 Du v2 Þn dt þ ðstuff in v1 ; v2 planeÞ
1 djdAj jdAj dt
We will also use d(/) where the Dirac delta function
where Du is the full derivative of u. Thus we obtain dðdAÞ ¼ ðDu ds1 Þ ds2 þ ds1 ðDu ds2 Þ dt
ð26Þ
So that using our above result yields
ð18Þ
For any tangent vector ds, we have dðdsÞ ¼ change in u over segment ds ¼ Du ds dt
dðdAÞ ¼ jdAjðr u n Du nÞ dt
ð14Þ
ð25Þ
where $ is the usual 3D divergence. Thus we get an expression only involving the normal of the surface, as desired
ð33Þ
½H ð/Þt ¼ H 0ð/Þ/t ¼ H 0 ð/Þðu r/Þ ¼ H 0 ð/Þjr/ju n This yields Z dV ¼ H 0 ð/Þjr/ju n dX dt X
ð34Þ
ð35Þ
j ¼ 0, then we must To preserve volume, if we have dV dt t¼0 have dtd ðdV Þ ¼ 0 dt Z d dV d H 0 ð/Þðui @ i /ÞdX ¼ dt dt dt X Z Z ¼ d0 ð/Þ/t ui @ i / dX H 0 ð/Þðui Þt @ i / dX X X |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} I II Z dð/Þui @ i /t dX ð36Þ X |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} III
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
1 dð/Þ u ru þ ð½pn rjn þ f þ FÞ r/ dX l X Z Z 1 dð/Þjr/jdX þ ðu ruÞ ndð/Þjr/jdX ¼ ½p l X ZX 1 þ ð37Þ ðrj f n F nÞdð/Þjr/jdX l Z X ðIIIÞ ¼ dð/Þui @ i /t dX Z ZX @ i ðdð/Þui /t Þ dX d0 ð/Þ@ i /ui /t dX ¼ X X Z dð/Þ@ i ui /t dX Z ZX ¼ d0 ð/Þ@ i /ui /t dX ðr uÞðu nÞdð/Þjr/jdX ðIIÞ ¼
Z
X
X
ð38Þ Therefore, Z Z d dV 1 dð/Þjr/jdX ðu ruÞ ¼ ½p dt dt X l X Z ndð/Þjr/jdX þ ðr uÞðu nÞdð/Þjr/jdX X Z 1 ðrj f n F nÞdð/Þjr/jdX ð39Þ X l From our volume preserving constraint, we can get [p]. R ½p ¼
X
ððr uÞðu nÞ ðu ruÞ n l1 ðrj f n F nÞÞdð/Þjr/jdX R 1 dð/Þjr/jdX X l
ð40Þ
2.2.3. Specification of l, dl and f dt To completely specify the equations of motion, l; dl (condt vective derivative of l) and f (additional momentum sources, including mass transfer effects) are required. No matter how these are specified, the motion will conserve volume and momentum, since the equations were derived from those principles. However, a physical bubble R motion should also conserve the mass of bubble (M ¼ surface ljdAj), unless we are modeling a process like evaporation or condensation, which can alter the mass. There are four cases for l and f. Each of these cases also has simple restriction to 2D. Case 1. l=constant on the surface and in time (e.g. l = 1, say)This does not conserve bubble mass, since mass will vary in proportion to the area, which varies in time. Thus, it does not correspond to a physical bubble motion. Since it is artificial, we may as well also take f = 0, for simplicity. This choice results in the simplest set of equation R ððn Du nÞu ðu ruÞÞ ndð/Þjr/jdX R ½p ¼ rj þ X dð/Þjr/jdX X ð41Þ
where j is the average curvature over the surface area.
527
Thus, the final equation of motion becomes ut þ u ru ¼ rðj jÞ uðr n Du nÞ Z 1 þ n ðn Du nÞðu nÞdð/Þjr/jdX A X Z 1 ð42Þ n ðu ruÞ ndð/Þjr/jdX A X Case 2. l = constant over the surface, but varies in time to M conserve mass(l ¼ AðtÞ and M = 1, say) This is the simplest motion that conserves mass, and corresponds to a real bubble with a fluid surface. So the mass is mobile enough to stay uniformly distributed. In this case there is mass transfer to a surface element during the motion. Since dtd ðljdAjÞ is not zero, and we need to postulate a form of the additional source of momentum f. We imagine that when mass transfers into/out of the element being considered, it takes its momentum with it. This suggests that the force FMT on the element ds, as described above, is proportional to the local rate of change of element mass d ðljdAjÞ. Also, since we are considering here an dt idealized mass transfer that is so rapid it keeps the mass density uniform on the surface, it in some sense samples and averages things over the whole surface, so that the mass undergoing transfer carries with it the average velocity of the entire bubble. This intuition leads one to postulate that the additional momentum transfer to the surface element ds due to the mass transfer is F MT ¼ u
d ðljdAjÞ dt
ð43Þ
where u is the average vector velocity of the bubble Z 1 ujdAj ð44Þ u¼ A MT and, since f ¼ FjdAj ,
1 d ðljdAjÞ jdAj dt 1 dl ¼ lu r u n Du n þ l dt
f ¼u
ð45Þ
In this case, we can again make the same simplifications as in Case 1 to get [p], Z 1 ½p ¼ rj þ 2 ðu uÞ nðr u n Du n A X Z 1 dl 1 dð/Þjr/jdX þ 2 ððu nÞðr uÞ þ l dt A X ðu ruÞ nÞdð/Þjr/jdX ð46Þ and the equation of motion is
528
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
lðut þ u ruÞ ¼ rðj jÞn þ lðu uÞ r u n Z 1 dl n Du n þ 2 ðu uÞ l dt A X 1 dl n r u n Du n þ l dt Z n dð/Þjr/jdX 2 ððu nÞðr uÞ A X ðu ruÞ nÞdð/Þjr/jdX ð47Þ Case 3. ljdAj is constant during the motion (i.e. d ðljdAjÞ ¼ 0).This also conserves mass; it corredt sponds to an elastic membrane, where the mass is not free to migrate. In this case, we envision no mass transfer effect, thus f = 0. We get, from d ðljdAjÞ ¼ 0, that dt 1 dl 1 djdAj ¼ ¼ ðr u n Du nÞ l dt jdAj dt
ð48Þ
Therefore, R ððr uÞðu nÞ ðu ruÞ n l1 rjÞdð/Þjr/jdX X R1 ½p ¼ dð/Þjr/jdX l ð49Þ
Case 4. l flows with some mass preserving velocity on the surface. This velocity would generally require additional equations to specify it. The local momentum source due to mass transfer would be f = $s(lul), the surface divergence of momentum flux associated with l, velocity ul.
dA ! dL jdAj ! jdLj A!L Note that dL = njdLj is an outward normal. 3.3. Initial data In order to have volume preserving motions, the initial velocity of the bubble must also be volume preserving. ¼ 0attime ¼ 0). From Eq. (34), we see this means (i.e. dV dt that the initial velocity u0 must satisfy Z 0¼ u0 dA ð50Þ surface
If this condition is violated, there will be very large acceleration produced when the time evolution starts. The easiest way to insure (50) is to either take u0 = 0 which corresponds to a bubble at rest initially, or u0 Æ n = 0 which corresponds to a bubble ‘‘rotating’’ initially (i.e., velocity tangent to surface everywhere). 3.4. General equation of bubble motion The general equation for bubble motion is du 1 djdAj 1 dl þ l ¼ ½pn rjn þ f lu dt jdAj dt l dt
ð51Þ
where 1 jdAj d ¼ r u n Du n jdAj dt
ð52Þ
and R
3. Summary of bubble motion equations We summarize our equations of bubble motion for various cases. All equations are valid in both 2D and 3D. 3.1. Notational convention We review the notations used below: dA
the vector area of surface element in the direction normal to the surface jdAj scalar area of the surface element n the unit outward normal to the surface j mean curvature r surface tension dQ ¼ Q þ t R u rQ for any quantity Q dt R Qave ¼ A1 surface QjdAj where A is the total area (A ¼ jdAj) 3.2. 2D versions of equation The 2D versions of all equations are obtained by replacing area of bubble surface by length of the bubble curve in all relevant places:
½p ¼
X
ððr uÞðu nÞ ðu ruÞ n l1 ðrj f n F nÞÞdð/Þjr/jdX R 1 dð/Þjr/jdX Xl
ð53Þ
In addition, we specify conservation of mass Z Mass : M ¼ ljdAj ¼ constant in time
ð54Þ
and no self-induced force Z 0 ¼ f jdAj
ð55Þ
3.5. Volume preserving acceleration motion Consider du ¼ ½pn rjn dt
ð56Þ
This preserves the bubble volume if Z 1 ½p ¼ rj þ ððr uÞðu nÞ ðu ruÞ nÞdð/Þjr/jdX A X ð57Þ
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
This model is volume preserving, but it does not conserve momentum or mass for the bubble, so it is not a physical bubble motion. However it is useful to compare it with volume preserving velocity case. 3.6. Volume and momentum preserving model
529
This preserves bubble volume, mass and momentum, so it is a physically realistic bubble motion. It corresponds to a bubble modeled as an elastic membrane – i.e. no migration of mass within the surface. Various interesting effects can be created by starting with non-uniform mass distributions.
This is the model described in Section 2.2.3, Case 1, with l=1
4. Numerical method
du ¼ ½pn rjn uðr u n Du nÞ dt
4.1. Description of the level set approach
where ½p ¼ rj þ
R X
ð58Þ
ððn Du nÞðu nÞ ðu ruÞ nÞdð/Þjr/jdX R dð/Þjr/jdX X ð59Þ
This model does not conserve bubble mass, but it does conserve volume and momentum.
We construct a level set function / such that bubble interface is the zero level set of /. We also initialize / to be the signed distance from the interface such that / is positive inside the bubble and negative outside the bubble. This is easy to do using the method in [6]. We represent the interface by surface of bubble ¼ fxj/ðxÞ ¼ 0g /ðxÞ > 0 inside bubble
ð68Þ
/ðxÞ < 0 outside bubble
3.7. Fluid bubble This is the model described in Section 2.2.3, Case 2, with bubble mass M = 1 du 1 dl l ¼ ½pn rjn þ f lu r u n Du n þ ð60Þ dt l dt 1 ð61Þ lðtÞ ¼ AðtÞ 1 dl f ¼ lu r u n Du n þ ð62Þ l dt Z 1 ½p ¼ rj þ 2 ðu uÞ n r u n Du n A X Z 1 dl 1 þ dð/Þjr/jdX þ 2 ððu nÞðr uÞ l dt A X ðu ruÞnÞdð/Þjr/jdX ð63Þ This preserves bubble volume, mass and momentum and so is physically realistic bubble motion. It corresponds to a bubble with an ideal fluid surface on which mass is very mobile and maintains a uniform mass distribution. 3.8. Elastic membrane bubble This is the model described in Section 2.2.3, Case 3 du ¼ ½pn rjn ð64Þ dt lt¼0 ¼ a positive valued function on surface ð65Þ 1 dl ¼ ðr u n Du nÞ ð66Þ l dt R ððr uÞðu nÞ ðu ruÞ n l1 rjÞdð/Þjr/jdX X R1 ½p ¼ dð/Þjr/jdX l l
ð67Þ
The idea of the level set method is to move / with the correct speed u at the front using the following differential equation: /t þ u r/ ¼ 0
ð69Þ
Next, we reinitialize, using the algorithm [6]. keeping / to be signed distance, at least near the front. Additionally, we save computational time performing these calculations only near the front. There are several localization algorithms available; we use the relatively simple algorithm developed in [7]. 4.2. Notations on using level set function We can rewrite the variables by using the level set function /. r/ jr/j r/ j ¼ r jr/j RR jdð/Þjr/jdx dy j¼ RR dð/Þjr/jdx dy 8 > <1 n¼
ð70Þ ð71Þ ð72Þ if / > a
0 if / < a px > :1 / 1 otherwise 1 þ þ sin a p a ( 2 1 if j/j < a 1 þ cos p/ a 2a
H a ð/Þ ¼
da ð/Þ ¼ Z Zsurface volume
0
QjdAj ¼
Z
otherwise
Qdð/Þjr/jdX ZX Q dV ¼ QH a ð/ÞdV X
ð73Þ
ð74Þ ð75Þ ð76Þ
530
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
Osher and Sethian [4]. Consider the following differential equation:
where a is the prescribed ‘‘thickness’’ of the interface (usually 3 grid points in our calculations) and Q is any quantity.
/t þ u r/ ¼ 0
4.3. Outline of the numerical method
For u a given function of space and time, we obtain higher order accuracy by using ENO type schemes both in time and in space. For the time discretization, we use Heun’s method [3]. Let
We can now summarize our algorithm. Step 1. Initialize / (x, t) such that / is a signed distance function to the front. Step 2. Solve the governing equation and get the velocity u and update the level function /. The method for updating the level function / was proposed by 1
1
0
1
0 t = 0.0125
1
0 t = 0.025
1
1
0
0
0
1
1
0 t = 0.05
1
1
0
1
0 t = 0.0625
1
0 t = 0.1
1
1
1
0
0
0
0 t = 0.025
1
0 t = 0.05
1
1
1
0
0
0
0 t = 0.075
1
1
1
0
0
0
0 t=0
1
0 t = 0.025
1
1
1
1
0
0
0
1
0 t = 0.075
1
0 t = 0.1
1
1
1
1
0
0
0
1
0 t = 0.1
1
0 t = 0.125
0 t = 0.15
1
0 t = 0.175
1
1
1
1
0
0
0
1
0 t=0
1
0 t = 0.0125
1
1
1
1
0
0
0
1
0 t = 0.0375
1
0 t = 0.05
1
1
1
1
1
1
0
0
0
0
0
0
1
1
0 t = 0.125
1
0 t = 0.2
1
Fig. 3. Area preserving curvature dependent velocity: spiral relaxes towards circle.
1
0 t = 0.15
0 t = 0.05
0
0 t = 0.0875
1
1
1
Fig. 1. Area preserving curvature dependent velocity: circle is an equilibrium solution.
0 t=0
Then we apply Heun’s method.
1
0
0 t = 0.075
ð78Þ
0
1
0 t = 0.0375
o / ¼ L½/ ot
1
0
0 t=0
ð77Þ
0 t = 0.175
1
0 t = 0.2
1
Fig. 2. Area preserving curvature dependent velocity: ellipse relaxes to circle.
0 t = 0.075
1
0 t = 0.0875
1
0 t = 0.025
1
0 t = 0.0625
1
0 t = 0.1
1
Fig. 4. Area preserving curvature dependent velocity: square relaxes to circle.
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
/mþ1 ¼ /m DtL½/m 1 1 Dt /mþ1 ¼ /m þ /mþ1 L½/mþ1 ð79Þ 2 2 2 For the space discretization, first consider the motion of curvature dependent velocity.We can rewrite Eq. (79) as /t ðj jÞn r/ ¼ 0
When we apply the numerical scheme to a Hamilton–Jacobi equation, care is required. For the term corresponding jjr/j, we use ENO scheme(we will explain it at step 3). All the derivatives of the term jj$/j are approximated by central differences. The term jj$/j depends on the 1 multiplication and/or division by ð/2x þ /2y Þ2 , which might be close to zero near some points. This implies that we should be consistent with an approximation to j and j$/j in the term jj$/j to avoid unnecessary errors. For the curva-
ð80Þ
Using the relation (72), (80) becomes /t þ jjr/j ¼ jjr/j
531
ð81Þ
1
1
1
1
1
1
0
0
0
0
0
0
0 t=0
1
0 t = 0.00625
1
0 t = 0.0125
1
0 t=0
1
0 t = 0.00625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.01875
1
0 t = 0.025
1
0 t = 0.03125
1
0 t = 0.01875
1
0 t = 0.025
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.05
1
Fig. 5. Area preserving curvature dependent velocity: starfish relaxes towards circle.
0 t = 0.0375
1
0 t = 0.04375
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.00625
1
0 t = 0.0125
1
0 t=0
1
0 t = 0.00625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.01875
1
0 t = 0.025
1
0 t = 0.03125
1
0 t = 0.01875
1
0 t = 0.025
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.05
1
Fig. 6. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
1
0 t = 0.03125
1
0 t = 0.05
1
Fig. 7. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
1
0 t=0
0 t = 0.0125
0 t = 0.0375
1
0 t = 0.04375
1
0 t = 0.0125
1
0 t = 0.03125
1
0 t = 0.05
1
Fig. 8. Area preserving curvature dependent velocity: merging case, relaxes towards circle.
532
M. Kang et al. / Computers & Fluids 37 (2008) 524–535 Table 1 Ratio of area for final area to initial area
Fig. 1
Fig. 3
Fig. 5
Mesh size
32 · 32
64 · 64
128 · 128
Initial area Final area Final/initial Initial area Ffinal area Final/initial Initial area Final area Final/initial
0.7873 0.7788 0.9892 0.9762 0.9343 0.9570
0.7861 0.7815 0.9942 0.9747 0.9580 0.9829 0.9261 0.9241 0.9979
0.7856 0.7832 0.9969 0.9740 0.9688 0.9947 0.9263 0.9249 0.9985
Table 2 Ratio of area for final area to initial area
Fig. 6
Fig. 7
Fig. 8
Mesh size
128 · 128
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Ffinal area Final/initial
1.1000 1.0810 0.9827 0.7224 0.6955 0.9627 0.8806 0.8928 0.9424
Fig. 9. Three dimensional volume preserving curvature dependent velocity: symmetric pinch off, relaxes towards two spheres. 1
1
1
0
0
0
0 t=0
1
0 t = 0.0125
1
1
1
1
0
0
0
0 t = 0.0375
1
0 t = 0.05
1
1
1
1
0
0
0
0 t = 0.075
1
0 t = 0.0875
1
0 t = 0.025
1
0 t = 0.0625
1
0 t = 0.1
1
Fig. 11. Curvature dependent acceleration (volume preserving): circle remains circle.
Fig. 10. Three dimensional volume preserving curvature dependent velocity: non-symmetric pinch off, relaxes towards two spheres.
ture dependent acceleration, first we solve the governing equation to get the velocity u. After that we solve the Eq. (77) using the second-order ENO scheme[5] for the space. For the x-direction in 2D(or 3D),
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
8 ui;j ðDx /i;j þ h2 minmodðDx Dxþ /i;j ; Dx Dx /i;j Þ > > > < if u > 0 i;j u r/ ¼ > ui;j ðDxþ /i;j h2 minmodðDxþ Dxþ /i;j ; Dxþ Dx /i;j Þ > > : otherwise minmodðx; yÞ ¼
y-direction and the z-direction for 3D problems also. Step 3. Construct a new distance function /.Consider the following equation: d t ¼ signðd 0 Þð1 jrdjÞ with d 0 ¼ /n
ð82Þ
1
0
1
0
0 t=0
1
1
1
1
0
0
0
0
0 t = 0.125
1
0 t = 0.25
1
0 t=0
1
0 t = 0.0625
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.375
1
0 t = 0.5
1
0 t = 0.625
1
0 t = 0.1875
1
0 t = 0.25
1
1
1
1
1
1
1
0
0
0
0
0
0
0 t = 0.75
1
0 t = 0.875
1
0 t=1
1
Fig. 12. Curvature dependent acceleration (volume preserving): oscillating ellipse.
0 t = 0.375
1
0 t = 0.4375
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.125
1
0 t = 0.25
0 t=0
1
1
1
1
1
0
0
0
0
0 t = 0.375
1
0 t = 0.5
1
0 t = 0.625
1
0 t = 0.0125
1 1
0
0 t = 0.0375
1
1
1
0 t = 0.05
1
1
1
1
1
1
0
0
0
0
0
0
1
0 t = 0.875
1
0 t=1
1
Fig. 13. Curvature dependent acceleration (volume preserving): oscillating square.
1
0 t = 0.3125
1
0 t = 0.5
1
0 t = 0.075
1
0 t = 0.025
1
0 t = 0.0625
1
0
1
0 t = 0.75
0 t = 0.125
Fig. 14. Curvature dependent acceleration (volume preserving): oscillating starfish.
1
0 t=0
ð84Þ
where sign(d0) gives the sign of d0. Given initial data, we solve the above equation until the solution reaches a steady state near the front. To eliminate the stiffness of sign function, we approximate sign(/) by
signðxÞ minðjxj; jyjÞ if x y > 0 ð83Þ 0 otherwise
and Dx ; Dxþ are just backward and forward divided differences. Similarly, we can apply this to the 1
533
0 t = 0.0875
1
0 t = 0.1
Fig. 15. Curvature dependent acceleration (volume and momentum preserving): circle remains circle.
534
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
/ signð/Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi / 2 þ 2
ð85Þ
where is very small number (i.e. = h). Following [3,6], we use the approximation: 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 > > maxððaþ Þ ; ðb Þ Þ þ maxððcþ Þ ; ðd Þ Þ > > > > > < if / > 0 0 jr/j ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > > maxðða Þ2 ; ðbþ Þ2 Þ þ maxððc Þ2 ; ðd þ Þ2 Þ > > > > : otherwise ð86Þ h minmodðDx Dxþ /i;j ; Dx Dx /i;j Þ 2 h b ¼ Dxþ / minmodðDxþ Dxþ /i;j ; Dxþ Dx /i;j Þ 2
a ¼ Dx / þ
1
1
1
0
0
0
0 t=0
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.5
1
1
1
1
0
0
0
0 t = 0.75
1
0 t = 0.875
1
0 t = 0.25
1
0 t = 0.625
1
0 t=1
1
Fig. 17. Curvature dependent acceleration (volume and momentum preserving): oscillating square.
a ¼ minða; 0Þ bþ ¼ maxðb; 0Þ b ¼ minðb; 0Þ
1
1
1
0
0
0
ð88Þ
We can define c, d using same approximation in the y-direction. For the criterion for steady state, we use Q¼
0 t = 0.125
1
ð87Þ
aþ ¼ maxða; 0Þ
1
n Rjd ni;j j
M
< ðDtÞh2
0 t=0
1
0 t = 0.0625
1
1
1
1
0
0
0
0 t = 0.125
1
0 t = 0.3125
1
0 t = 0.5
1
ð89Þ jd ni;j j
where M = number of grid point where < að¼ chÞ, for some constant c (we use c = 1.5). This
1
1
1
0
0
0
0 t=0
1
0 t = 0.125
1
0 t = 0.25
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.5
1
1
1
0
0
0
0 t = 0.75
1
0 t = 0.875
1
1
0 t = 0.25
1
1
1
1
0
0
0
0 t = 0.375
1
0 t = 0.4375
1
Fig. 18. Curvature dependent acceleration (volume and momentum preserving): oscillating starfish.
is a very local criterion and convergence is very rapid.After solving this equation we let 0 t = 0.625
1
1
0 t = 0.1875
1
/nnew ¼ steady state solution of ð84Þ:
ð90Þ
Step 4. We have now advanced one time step. Go to step 2 and repeat.
0 t=1
1
Fig. 16. Curvature dependent acceleration (volume and momentum preserving): oscillating ellipse.
5. Numerical results in 2D and 3D Our experiments simulate the motions of soap bubbles. In the case of curvature dependent velocity and accelera-
M. Kang et al. / Computers & Fluids 37 (2008) 524–535
tion, we simulate the merging and breaking of the interface in 2D and 3D. 5.1. Results using curvature dependent velocity We simulate several different shapes in the case of curvature dependent velocity. In Figs. 1–5, we show the continuous evolution of non-intersecting curves collapsing smoothly to a circle. In Figs. 6–8, we show the merging of two bubbles in 2D. After merging they converge to a steady state which is, of course, a circle. Table 3 Ratio of area for final area to initial area
Fig. 11
Fig. 12
Fig. 13
Fig. 14
Mesh size
32 · 32
64 · 64
128 · 128
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial
0.7873 0.7806 0.9915 0.8792 0.7939 0.9030 1.0000 0.9856 0.9856 0.8873 0.7964 0.8976
0.7861 0.7831 0.9962 0.8790 0.8413 0.9570 1.0000 1.0240 1.0240 0.9201 0.8470 0.9205
0.7856 0.7841 0.9980 0.8791 0.8635 0.9823 1.0000 1.0320 1.0320 0.9267 0.8882 0.9582
Table 4 Ratio of area for final area to initial area
Fig. 15
Fig. 16
Fig. 17
Fig. 18
Mesh size
32 · 32
64 · 64
Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial Initial area Final area Final/initial
0.7873 0.7806 0.9915 0.8792 0.7927 0.9017 1.0000 0.9810 0.9810 0.8873 0.8105 0.9135
0.7861 0.7831 0.9962 0.8790 0.8414 0.9572 1.0000 1.0130 1.0130 0.9201 0.8744 0.9503
535
In Figs. 9 and 10, we use a symmetric and non-symmetric dumbbell as initial data. In 2D, a dumbbell under our motion collapses to a circle. However, in 3D, a dumbbell breaks into two part and each part changes smoothly to a sphere. Both 2D and 3D motions preserve the area and the volume, respectively (see Tables 1 and 2). 5.2. Results using curvature dependent acceleration In Figs. 11–14, we simulate the motion of bubble which was previously described in Section 2.2.3, Case 1. The result shows that motion has periodic oscillations as expected. We next use the motion described in Section 2.2.3, Case 2 (see Figs. 15–18). Both cases preserve area well (see Tables 3 and 4). Acknowledgements M. Kang was supported by Korea Research Foundation Grant (R08-2004-000-10233-0). M. Kang also holds joint appointment in the Research Institute of Mathematics, Seoul National University. S. Osher was supported by ONR Grant N00014-03-1-0071. References [1] Chang YC, Hou TY, Merriman B, Osher S. A level set formulation of Eulerian interface capturing methods for incompressible fluid flows. J Comput Phys 1996;124:449–64. [2] Hou TY, Lowengrub J, Shelley M. Removing the stiffness interfacial flows with surface tension. J Comput Phys 1994;114:312–38. [3] Osher S, Shu C-W. High-order essentially nonoscillatory schemes for Hamilton–Jacobi equation. SINUM 1991;28:907–22. [4] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed : algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79:12–49. [5] Shu C-W, Osher S. Efficient implementation of essentially nonoscillatory shock capturing schemes II. J Comput Phys 1989;83:32– 78. [6] Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. J Comput Phys 1994;114:146–59. [7] Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE based fast local level set method. J Comput Phys 1999;55:410–38.