有介面活性劑之曲率運動的數值方法

國立交通大學

應用數學系數學建模與科學計算碩士班

碩士論文

有介面活性劑之曲率運動的數值方法

Numerical methods for motion by curvature with

surfactant

研 究 生：許哲維

哲維

指導教授：賴明治 教授

有介面活性劑之曲率運動的數值方法

學生：許哲維 指導教授：賴明治

國立交通大學

應用數學系數學建模與科學計算碩士班

摘要

在這篇論文裡，我們使用兩種數值方法去處理曲率運動的問題。

第一種方法稱為介面追蹤法(front-tracking method)，主要是用有限差

分法追蹤曲率運動後的曲線。為了瞭解非線性表面張力的影響，我們

將介面活性劑加至曲線上，讓介面活性劑能夠在曲線上自由的擴散。

我們所介紹的數值方法能夠保證隨著時間的流逝介面活性劑之濃度

的總和不變。另一個方法稱為多重等位函數法(multiple level set

method)，主要用此方法處理圖形結構變化的問題。圖形上的每條曲

線的運動是由曲率運動所決定。在此問題中我們暫時不添加介面活性

劑。

Numerical Methods for Motion by

Curvature with Surfactant

Student : Che-Wei Hsu Advisor : Ming-Chih Lai

Department of Applied Mathematics

Institute of Mathematical Modeling and Scientific Computing

National Chiao Tung University

Abstract

In this paper, we use two numerical methods to model a grain

boundary evolution for motion by mean curvature. In the first part,

we present a finite difference method to track a network of curves

whose motion is determined by curvature. To study the effect of

inhomogeneous surface tension on the evolution of the network of

curves, we include surfactant which can diffuse along the curves.

Our numerical method is based on a direct discretization of the

governing equations which conserves the total surfactant mass in

the curve network. In the second part, we present a multiple level

set method to track the grain topological change. The motion of

grain boundary is considered by mean curvature. In this part, there

is no surfactant on each grain boundary. The governing equations

consist of a level set equation for the grain boundary motion and

reinitialization equation for reconstructing signed distance function.

Thus, one important step is coupling of multiple level set functions.

Numerical examples shows the topological change evolution of the

grain structure.

誌 謝

在建模所兩年的學習過程，隨著論文的付梓，即將劃上句點，這

段時間以來的點點滴滴，有回憶，有不捨；回憶之情將在我的懷中日

漸晶瑩光耀，不捨之心將使我的人生成就勇氣。

本論文能順利完成，首先要感謝我的指導老師—賴明治教授，對

於研究的方向、觀念的啟迪、程式的指導、參考資料的提供與求學的

態度逐一斧正與細細關懷，讓我漸漸熟悉數值偏微分方程與科學計算

等領域，於此獻上最深的敬意與謝意。論文口試期間，承蒙口試委員

楊肅煜教授與吳金典助理教授的鼓勵與疏漏處之指正，使得本論文更

臻完備，在此謹深致謝忱。

在研究所修業期間，科學計算實驗室的學長們給我許多的幫助。

感謝陳冠羽學長對於數值方法的指點，並且提供許多參考資料給我，

讓我在研究的過程中能夠持續獲得進展。感謝曾昱豪學長、黃仲尹學

長以及曾鈺傑學長不僅在研究上提供我保貴的意見，更是難得的球友，

讓我在球場上能盡情的釋放壓力，處於最佳的備戰狀態。感謝建模所

的好夥伴們─昆霖、仁洲、裕昇、昱丞以及振庭，有你們的陪伴，研

究也能伴隨著笑聲，讓我的心情輕鬆愉快。

最後，特將本文獻給我最敬愛的母親，感謝您無怨無悔的養育與

無時無刻的關懷照顧，還有父親及弟弟哲瑞在經濟上與精神上的支持，

讓我能專注於課業研究中，願以此與家人共享。

目 次

2 The governing equations for front-tracking method 2

3 Numerical method 4

4 Numerical results 8

4.1 Three-curve network . . . 8

4.2 von Neumann law . . . 9

5 The governing equations for multiple level set method 14

6 Numerical method 16

7 Numerical results 18

7.1 Types of grain topological change . . . 18

7.2 von Neumann law . . . 20

7.3 Imposition of periodic boundary conditions . . . 21

1 Introduction

In this paper, we generalize the grain growth problem discussed in [1] by including the eect of an inhomogeneous surface tension. For practical problems, it is dicult to maintain constant surface tension as insoluble surface active agents (surfactant) are common and their presence could signicantly aect the value of the surface tension, therefore the dynamics of interface motion [6]. To account for the eect of the surface tension on the interfacial dynamics of a complex network of interfaces, we consider a network of curves in two dimensional setting and assume that there is surfactant distributed along the curve and the surface tension varies according to the surfactant concentration.

In recent years the processes of grain boundary evolution and grain growth have been studied extensively through computer simulations. Direct computer simulations of the grain structure and topological changes can be further classied into proba-bilistic and deterministic approaches. Probaproba-bilistic models include the Potts model [7],[8] and kinetic lattice Monte Carlo methods [9]. In deterministic models a precise description of the motion of the grain boundaries is needed. It is assumed that the grain boundary evolution velocity is proportional to the driving force. Front-tracking method [1],[10] use this assumption. In front-tracking method, the grain boundaries are discretized and tracked explicitly by piecewise curved segments. However, level set method use an implicit representation of the geometry.

The level set method was rst proposed in 1980 [5] for front propagation with cur-vature dependent speed. Since than, it has been extended to numerous applications with moving interfaces in uid mechanics, computer animation, image processing, among others, see [4]. In the level set method the interface is represented as the zero contour of a level set function which satises an evolution equation. The physics which governs the motion of the interface is included naturally in this model.

The rest of the paper is organized as follows. In Section 2, we present the gov-erning equations for front-tracking method which includes one parabolic equation for the curve motion coupled with a convection-diusion equation for the surfactant con-centration along each curve. The numerical method is described in Section 3 which includes an algorithm of solving the parabolic equation and a conservative scheme for the surfactant concentration equation. The eects of inhomogeneous surface tension on the motion of network and the cell evolution of the von Neumann law are inves-tigated numerically in Section 4. In Section 5, we present the governing equations for multiple level set method which includes a level set equation for the grain bound-ary motion and reinitialization equation for reconstructing signed distance function. The numerical method is described in Section 6. The topological change evolution of the grain structure is investigated numerically in Section 7. Some conclusions and remarks are given in Section 8.

2 The governing equations for front-tracking method

As in [1], we consider the situation with three phase boundaries in a unit circular domain described by parametric curves Xi_{(s; t), s 2 [0; 1], for i = 1; 2; 3 as shown in}

Fig. 1. The mean curvature motion of equations are dened as Xi t = i X i ss jXi sj2; (2.1)

where i_{(s; t) is the surface tension along the curve X}i _{and is determined by}

surfac-tant concentration i_{(s; t). In this paper we use the simplied nonlinear Langmuir}

equation of state [17]

i _{= }

c(1 + ln(1 i i)); (2.2)

where c is the surface tension of a clean interface and i satisfying 0  i < 1,

8i are dimensionless numbers that measure the sensitivity of surface tension to the surfactant concentration. Along the curve, the surfactant concentration is governed by the transport-diusion equation [2, 19]

i

t + (rs
Ui) i =  r2s i; (2.3)

where rs is surface gradient operator and r2s = rs
rs is the surface Laplacian.

Ui _{is the velocity of the curve and  is the diusion coecient of the surfactant}

concentration. In this paper, the diusion coecient of the surfactant concentration is constant for all curves. Throughout this paper, we dene

i_{(s; t) =} Xsi

jXi

sj (2.4)

as the unit tangent vector of the curve i. For simplicity, the above surface divergence and surface Laplacian can be written explicitly as

rs
Ui = @U i @i
i =  @Ui @s
i   @X_@si ; (2.5) r2 s i = _@s@  @ i @s  @X_@si @X_@si : (2.6) At s = 0, we assume that each curve meets the domain boundary with a normal angle. To be more precise, let b() be the given parametric representation of the unit circular domain with  2 [0; 2]. Then the conditions at the domain boundary are given by

Xi_{(0; t) = b(}i_{); (2.7)}

for all i such that

i_{(0; t)
b}0₍i_{) = 0: (2.8)}

We also assume that there is no surfactant ux across the domain boundary rs i(0; t)
i(0; t) = @

i_{(0; t)}

@s = 0: (2.9)

At triple junction (s = 1), three curves must meet at triple junction point, that is

X1_{(1; t) = X}2_{(1; t) = X}3_{(1; t): (2.10)}

In order to balance the surface tension at the triple junction, the condition is given by

3

X

i=1

i_{(1; t)}i_{(1; t) = 0: (2.11)}

The above condition comes from the Young-Laplace equation. Moreover, we also need to impose similar conditions for surfactant concentration at the triple junction. Firstly, the following condition hold the continuity of the surfactant concentration

1_{(1; t) =} 2_{(1; t) =} 3_{(1; t): (2.12)}

In order to balance the tangential ux for surfactant concentration at the triple junc-tion, the condition is given by

3 X i=1 rs i(1; t)
i(1; t) = 3 X i=1 @ i_{(1; t)} @s  @X_@si(1; t) = 0: (2.13) When i_{in Eq. (2.2) are identical for each curve, the surface tensions take same value}

the condition (2.11) becomes P3_i=1i_{(1; t) = 0, that is, the angles between these}

curves at the triple junction are 120 _{which results in the famous law of Plateau.}

Finally, we want to check that the total surfactant mass along the three curves is conserved. We rewrite the surfactant concentration equation in a form as

@ i @t  @X_@si +_@t@ @X_@si i _{= } @ @s  @ i @s  @X_@si; (2.14) by multiplying Eq. (2.3) with stretching factor j@Xi_{=@sj and using Eq. (2.5,2.6) and}

the following equation [2] (rs
Ui)  @X_@si = @U_@si
i ₌@Ui @s
@Xi @s   @X_@si =  @ @t  @Xi @s 
@X_@si @X_@si = _@t@ @X_@si: (2.15) From the domain boundary condition and the junction boundary condition, the total surfactant mass along the three curves is written in form as

3 X i=1 Z ₁ 0  @ i @t  @X_@si +_@t@ @X_@si i_{ds = 0: (2.16)}

By taking the time derivative outside the integral, we can get that the total surfactant mass along the three curves is conserved, i.e.,

d d t 3 X i=1 Z ₁ 0 i_{(s; t)}@Xi @s  ds ! = 0: (2.17)

3 Numerical method

For each parametric curve i, we set up a mesh sk = (k 1=2)
s, k = 1; 2; : : : m

where
s = 1=m, and use Xi

k = Xi(sk; n
t) as Lagrangian markers to represent the

curve at time t = n
t. The surfactant concentrations and surface tensions on each

curve are also dened on these Lagrangian markers and denoted by i

k = i(sk; n
t)

and i

k = i(sk; n
t), respectively. In order to avoid confusing reader, we use the

variables with tilde as the values at the next time step; that is, ~Xi

k = Xi(st; (n+1)
t)

and ~i

k= i(sk; (n + 1)
t). The whole procedures of the numerical time integration

are as follows.

Step 1. Compute the surface tension on each curve i i

Figure 2: The implementation of the domain boundary condition.

where c is the surface tension of a clean interface and i satisfying 0  i < 1 for

all i are dimensionless numbers that measure the sensitivity of surface tension to the surfactant concentration.

Step 2. Solve the equation of motion (2.1) by a explicit scheme as in [1] ~ Xi k Xki
t = ki(X i k+1 2Xki + Xk 1i )=
s2 j(Xi k+1 Xk 1i )=(2
s)j2 : (3.2)

As the above discretization, we use the central dierence schemes to approximate the rst and second derivatives. Note that the interior points k = 1; 2; : : : m are computed by the above discretization only. Next we provide the details on how to

nd the boundary points ~Xi

0 and ~Xm+1i at the domain boundary s = 0 and the

triple-junction boundary s = 1, respectively.

(a) At the domain boundary, we discretize Eq. (2.7) by central dierence approx-imation as

~ Xi

1+ ~X0i

2 = b(i): (3.3)

If the point b(i_{) is known, the domain boundary point can be extrapolated by}

~ Xi

0 = 2b(i) X~1i: (3.4)

As in Fig. 2, the point b(i_{) is determined by discretizing Eq. (2.8) with central}

dierence approximation to obtain ~ Xi

1 b(i)

s
b0(i) = 0: (3.5)

Since the domain is the unit circle and ~Xi

1 has polar coordinates (r; ), we obtain

b(i_{) = (cos ; sin ).}

(b) At the triple-junction boundary, we discretize Eq. (2.10) as ~ X1 m+1+ ~Xm1 2 = ~ X2 m+1+ ~Xm2 2 = ~ X3 m+1+ ~Xm3 2 = ~Xj: (3.6)

Figure 3: The implementation of the triple-junction boundary condition. If the triple-junction point ~Xj is known, the points ~Xm+1i are determined by

~ Xi

m+1 = 2 ~Xj X~mi : (3.7)

The detail of location of the markers at the triple-junction can be found in Fig. 3. We discretize Eq. (2.11) by forward dierence approximation as

1 m ~ Xj X~m1 j ~Xj X~m1j + 2 m ~ Xj X~m2 j ~Xj X~m2j + 3 m ~ Xj X~m3 j ~Xj X~m3j = 0; (3.8)

and solve these equations. Then we can obtain the triple-junction point ~Xj.

Step 3. Update surfactant concentration ~i

k as follows. For clarity, we denote

the discrete stretching factor j@Xi_{=@sj by}

jDsXkij =  Xk+1i Xk 1i 2
s  : (3.9)

By using the above approximation, we obtain

jDsX_k+1=2i j=   Xik+3=2 X i k 1=2 2
s   =   (Xik+3=2+X i k+1=2)=2 (Xk+1=2+Xk 1=2)=2
s   =   Xik+1 X i k
s    and jDsX_{k 1=2}i j=   Xik+1=2 X i k 3=2 2
s   =   (Xik+1=2+X i k 1=2)=2 (Xk 1=2+Xk 3=2)=2
s   =   Xik X i k 1
s   :

Next, we can discretize Eq. (2.14) by an explicit and symmetric schemes as ~i k ki
t jDsX~kij + jDsXkij 2 + jDsX~kij jDsXkij
t ~i k+ ki 2 =  1
s ( i k+1 ik)=
s jDsX_k+1=2i j ( i k ik 1)=
s jDsX_{k 1=2}i j ! ; (3.10)

on the interior points k = 1; 2; : : : m. For the values on the boundary points ~i o

and ~i

m+1, we use the following conditions at the domain boundary s = 0 and the

triple-junction boundary s = 1, respectively.

(a) At the domain boundary s = 0, we use central dierence approximation to discretize Eq. (2.9) and obtain ~i

0 = ~1i.

(b) At the triple-junction boundary s = 1, we approximate Eqs. (2.12) by ~1 m+1+ ~m1 2 = ~2 m+1+ ~m2 2 = ~3 m+1+ ~m3 2 = ~j; (3.11)

where ~j represents the surfactant concentration at the triple-junction. Moreover, we

discretize Eq. (2.13) by nite dierence scheme as (~1 m+1 + ~m1)=
s jDsX~_m+1=21 j + (~2 m+1+ ~m2)=
s jDsX~_m+1=22 j + (~3 m+1+ ~m3)=
s jDsX~_m+1=23 j = 0: (3.12)

By substituting Eq. (3.11) into Eq. (3.12), we solve a single equation to obtain the surfactant concentration at the triple junction ~j. Once ~j is found, ~ji can be

obtained from Eq. (3.11).

Note that, by taking the summation of both sides of Eq. (3.10) and using the numerical boundary condition (3.12) at the triple-junction, one can verify that

3 X i=1 m X k=1 ~i kjDsX~kij
s = 3 X i=1 m X k=1 i kjDsXkij
s: (3.13)

This is the discrete version of the conservation of total surfactant mass along all the curves, corresponding to the mid-point rule discretization for the integral in Eq. (2.17). It is interesting to note that the numerical scheme (3.2) for Eq. (2.1) is independent of the mesh width
s. Since the scheme is explicit, the time step size must be chosen similarly in [1] by
t = 1₈min j;k jXi k+1 Xk 1i j2 i k : (3.14)

Under this constraint, the time step becomes smaller if the length of any curve short-ens in which the marker spacing becomes smaller. One way to maintain the marker resolution is to delete the markers in an appropriate way so that the time step size can be maintained. On the other hand, if the curve stretches and the marker spacing is too coarse, then we need to add more markers along the curve. One important thing during the marker redistribution process is to keep the mass conservation of the surfactant. This can be done in a local way. For instance, in the segment of adding more marker points, we simply distribute the surfactant into those points uniformly. On the other hand, in the segment of removing marker points, we add those surfac-tant to the new combining segment. Thus, the overall surfacsurfac-tant mass is conserved exactly.

4 Numerical results

As in [1], we consider the network of curves inside the unit circle such that b() = (cos ; sin ), where  2 [0; 2]. For comparison purpose, we will present results for the cases without surfactant (clean) and with surfactant (contaminated). In Eq. (2.2), i _{= 0 implies no surfactant exists on the curves, in which we do not need}

to solve the surfactant equation (2.3). Thus, the clean curves have uniform surface

tension  = c = 1. In the following case, we choose the mesh width
s = 1=16 on

each curve and the time step
t = 0:0001.

4.1 Three-curve network

In the rst example, we consider the evolution of three curves in the unit circle with initial conguration as in [1]

X1_{(s; 0) = (1 s)( 1=2;} p_{3=2); (4.1)}

X2_{(s; 0) = (1 s)( 1=2;}p_{3=2); (4.2)}

X3_{(s; 0) = (1 s; sin}2_{(s)=4): (4.3)}

The above initial conguration is shown in Fig. 1. In order to check the eect the surfactant on the curves, we compare the cases with (i _{= 0:25; i = 1; 2; 3) and}

with-out surfactant (i_{; i = 1; 2; 3). For the case with surfactant, the diusion coecient}

Figure 4: The time evolution of curve networks: solid lines for i _{= 0 and dotted}

Figure 5: Distribution of the surfactant concentration on curve 1 (left), curve 2 (center) and curve 3 (right).

is chosen as  = 0:1. The initial surfactant concentration is uniformly distributed only along the curve 3 so that

1_{(s; 0) =} 2_{(s; 0) = 0;} 3_{(s; 0) = 1: (4.4)}

Fig. 4 shows the time evolution of these three curve networks. We denote the clean curve network (without surfactant) by solid line, while the contaminated curve net-work (with surfactant) by dotted line. As demonstrated in [1], the curve 3 will be attered because of the motion of mean curvature. Thus, the area between curve 2 and curve 3 decrease as time evolves. Since the existence of the surfactant along the curve 3 reduces the surface tension, the curve motion becomes slower than the clean curve motion. Fig. 5 shows the distribution of the surfactant concentration on each curve. In Fig. 5, x-axis and y-axis represents the length of the curve and the surfactant concentration on each curve. Due to the eect of diusion, the surfactant on curves 1 and 2 are no longer zero.

4.2 Von Neumann law

In 1952, von Neumann [14] showed that the rate of the change of the area of a given bubble (a curve polygon) in two-dimensional dry foam is independent of bubble size and solely dependent on the number of edges. The original derivation is based on the rate of gas diuses through a permeable wall. In our case, the rate of area

change of the domain is given by dA dt = X i Z Ui
_ni_{ds =}Z P ds; (4.5)

Figure 6: The time evolution of six nodes cell. Solid line: without surfactant; Dotted line: with surfactant. 1_{(s; 0) =} 2_{(s; 0) = 0;} 3_{(s; 0) = 1 and rest of} i_{(s; 0) = 0;}

1 _{= 0:75;} 2 _{= 0:5;} 3 _{= 0:25 and the rest of} i _{= 0:25.}

Figure 7: The time evolution of ve nodes cell. Solid line: without surfactant; Dotted line: with surfactant. 1_{(s; 0) =} 2_{(s; 0) = 0;} 3_{(s; 0) = 1 and rest of} i_{(s; 0) = 0;}

Figure 8: The time evolution of seven nodes cell. Solid line: without surfactant;

Dotted line: with surfactant. 1_{(s; 0) =} 2_{(s; 0) = 0;} 3_{(s; 0) = 1 and rest of}

i_{(s; 0) = 0;} 1 _{= 0:75;} 2 _{= 0:5;} 3 _{= 0:25 and the rest of} i _{= 0:25.}

Figure 9: The cell area versus time, for the cases with and without surfactant. With-out surfactant (solid lines), the area evolves as predicted by the von Neumann law while for the case with surfactant (dotted lines), the rate of growth is increased for n = 7 while the rate of decay for n = 5 is decreased. The area for n = 6 with surfactant also increases.

where Ui
_ni _{represents the normal velocity of the curve i per unit length, as discussed}

earlier. Note that, the above integral of mean curvature is over all the curves that enclose the area. When  is a constant, the enclosed area is polyhedral-shaped with

arcs edges, the above integral can be simplied as dA dt =  n X i=1 i 2 ! = ₃ (n 6); (4.6)

where i = =3 is the exterior (turning) angle at the vertex, and n is the number of

edges. A similar derivative can be found in [12]-[14]. When surface tension varies,

however, Eq. (4.6) is no longer valid. For example, i does not always take the

value of =3. It will be interesting to examine how the area changes under a similar setup. We start with a single n-vertices inner cell with circular arcs and connect those vertices with n straight lines to the domain boundary. In particular, the cell boundary is a circle of radius 0:5. Note that, the number of lines is the same as the number of vertices on the inner cell. Along each line or circular arc, we lay out a parametrization on those curves. So a cell network with n vertices should have 2n curves and n triple-junction.

From the von Neumann law, one can see that the cell with more than six edges will grow while the one less than six edges will shrink when the surface tension is constant. More specically, the area of a cell with n = 6 remains unchanged while the cells with n = 5 and n = 7 should have the same growth (decay) rate. Fig. 6 shows the time evolution of a six vertices cell with (dotted line) and without surfactant (solid line). One can see that in the absence of surfactant, the cell area does not change. However, with surfactant, the system behaves dierently. For illustrative

purposes, we add surfactant only along one line segment 3_{(s; 0) = 1 initially, and}

choose 1 _{= 0:75;} 2 _{= 0:5 and} 3 _{= 0:25 in Langmuir equation along the rightmost}

three curves network (same set up as in Fig. 1) and keep other i _{= 0:25, then the}

symmetry is broken due to unbalanced surface tension, as shown in Fig. 6.

Figs. 7 and 8 show the time evolution of ve and seven vertices cell with (dotted line) or without surfactant (solid line), respectively. These gures have same set up as in Fig. 6. It is interesting to see that the unbalanced surface tension slows down the decay rate for the ve vertices cell area and speeds up the area growth for seven vertices cell. This is due to the fact with dierent i _{and surfactant concentration,}

thus, dierent surface tension along the rightmost three curves network. Fig. 9 shows the plot of cell area versus time for the cell with and without surfactant for dierent nodes n = 5; 6; and 7. One can see the result conrms the prediction by the von Neumann law when there is no surfactant.

As our nal example, we present the results when surfactant is added initially

i_{(s; 0) = 1 to all the outside line segments connecting the center network (cell) to}

the boundary. The inner cell boundaries are thus kept clean (without surfactant) initially. In this case, the initial exterior turning angles of the center network are all less than =3, which reduces the value of P7_i=1i 2 in the von Neumann law.

surfactant. Furthermore, the cell area should increase again as the surfactant diuses into the center network, as shown in Figs. 10 and 11.

Figure 10: The time evolution of seven nodes cell. Solid line: without surfactant; Dotted line: with surfactant. i _{= 0:25.}

Figure 11: The cell area versus time, for the cases with and without surfactant for n = 7. Without surfactant (solid lines), the area evolves as predicted by the von Neumann law while for the case with surfactant (dotted lines), the area decreases in the beginning but increases later on as the surfactant diuses into the center network.

5 The governing equations for multiple level set

method

We review the level set method as in [15]. Consider a moving interface X(t)

bounding a open region (t) in Rn _{of codimension one by a Lipschitz continuous}

function (x; t) which has the following properties 8 < : (x; t) < 0 if x 2 X(t) (x; t) = 0 if x 2 X(t) (x; t) > 0 if x 2 X(t)+ ; (5.1)

where X(t) and X(t)+ _{represents inside and outside of the interface respectively.}

This means that the interface X(t) is represented as the zero level set of function (x; t). We use the convection equation to dene the motion of the interface where (x; t) = 0 in a form as

t+ v
r = 0; (5.2)

where v is the interface velocity. We consider the motion by curvature where the interface moves in the normal direction with a velocity proportional to its curvature; i.e., v = bn, where b =  > 0 is a constant,  is the grain boundary mobility,  is the surface tension,  is the curvature and n is unit normal vector. Therefore, the interface moves in the direction of concavity, so that circles shrink to a single point and disappear. By plugging the interface velocity into the level set equation (5.2), we have t bn
r = 0: (5.3) Furthermore, since  = r
n = r
 r jrj  (5.4) and n
r = r jrj
r = jrj2 jrj = jrj; (5.5)

the level set equation for grain boundary motion induced by curvature is rewritten as

t b r
 r jrj  jrj = 0; (5.6) where jrj =p2

x+ 2y: In this paper, the Neumann boundary condition is used as

@

@n = 0 on @; (5.7)

where n is the unit normal vector on the domain boundary @. The initial conditions of (x; t) are often dened to be the signed distance function to the interface as

(x; 0) = 8 < : d (x; X(0)) if x 2 X(0) 0 if x 2 X(0) +d (x; X(0)) if x 2 X(0)+ ; (5.8)

Figure 12: The interface is constructed by four level set functions. where the distance function d (x; X) is dened as

d (x; X) = min(jx xIj) 8 xI 2 X: (5.9)

In addition, the signed distance function satises the following property

jrj = 1: (5.10)

For example, we use a signed distance function (x; y) = px2 _{+ y}2 _{1 to represent}

implicitly the unit circle. By computing x and y, we can get that

jrj =q2 x+ 2y = 2 4 _p x x2_{+ y}2 !₂ + p y x2 _{+ y}2 !₂3 5 1 2 = 1: (5.11)

In fact, the initial grain structure is described by multiple level set functions. For example, the initial grain structure as shown in Fig. 12 contains four level set func-tions. Each level set function i _{represents a grain. The grains change their topology}

by Eq. (5.6) independently. Therefore, in order to couple these grains, we do some corrections [15] for level set functions as

^i ₌ 1 2  i _min j6=i  j  ; (5.12)

where ^i _{is the corrector of}i_{. As we know, a number of simplications can be made}

when  is a signed distance function. For this reason, a reinitialization procedure is needed to reset the level set function to be a signed distance function to the interface. As in [4], we use the following reinitialization equation for the correction of  for all i at time t

@

@t = S(^)(1 jrj); (5.13)

where the sign function S(^) is given by S(^) = 8 > < > : 1 if ^ > 0 0 if ^ = 0 1 if ^ < 0 : (5.14)

6 Numerical method

Assume that the initial grain structure is composed of multiple signed distance functions i_{, i = 1; 2; : : : p. For each signed distance function i, we use the domain of}

two dimension D = [; ][ ; ], and set up the mesh of two dimension xk= +k
x,

k = 0; 1; : : : m and yh = + h
y, h = 0; 1; : : : m, where
x = ( )=m and

y = ( )=m. Thus, we use i

k;h = i((xk; yh); n
t) and ik;h = i((xk; yh); n
t)

to represent the grid point x = (xk; yh) on the domain D and the curvature at the

time t = n
t, respectively. In order to avoid confusing reader, we use the variables with tilde as the values at the next time step; that is, ~i

k;h= i((xk; yh); (n + 1)
t).

The whole procedures of the numerical time integration are as follows. Step 1. Compute the curvature of each grain boundary in a form as

i _{= r}
 ri

jri_j



= xxi (yi)2_[(i2 xiyixyi + yyi (xi)2

x)2+ (yi)2]3=2 : (6.1)

In this work, we approximate the rst and second partial derivative for spacial (i x,

i

y, xxi , yyi and xyi ) by central dierence approximation as

Dxk;hi = (k+1;hi k 1;hi )=(2
x); (6.2)

Dyk;hi = (k;h+1i k;h 1i )=(2
y); (6.3)

Dxxk;hi = (k+1;hi 2k;hi + k 1;hi )=(
x2); (6.4)

Dyyk;hi = (k;h+1i 2k;hi + k;h 1i )=(
y2); (6.5)

Dxyk;hi = (k+1;h+1i k 1:h+1i k+1;h 1i + k 1;h 1i )=(4
x
y): (6.6)

Next, the discretization of the curvature of each grain boundary at the point (xk; yh)

can be written by i

k;h=

Dxxk;hi (Dyk;hi )2 2Dxk;hi Dyk;hi Dxyk;hi + Dyyk;hi (Dxk;hi )2

[(Dxk;hi )2+ (Dyk;hi )2]3=2 : (6.7)

Step 2. Solve the level set equation (5.6) by explicit scheme in a form as ~i

k;h k;hi

t bk;hi [(Dxk;hi )2+ (Dyk;hi )2]

1

2 _{= 0; (6.8)}

where b is a constant. For the Neumann boundary conditions, we discretize Eq. (5.7) by central dierence and obtain the following equations

~i m;h= ~m 1;hi ; 8 h = 1; 2; : : : m 1; ~i 0;h = ~1;hi ; 8 h = 1; 2; : : : m 1; ~i k;0= ~k;1i ; 8 k = 1; 2; : : : m 1; ~i k;m= ~k;m 1i ; 8 k = 1; 2; : : : m 1; (6.9)

Figure 13: The deformation of each grain without coupling. and

~i

0;0 = ~1;1i ; ~m;0i = ~m 1;1i ; ~m;mi = ~m 1;m 1i ; ~0;mi = ~1;m 1i : (6.10)

Note that, there are multiple signed distance function in our initial grain structure. Thus, the deformation of each grain is independent if we do not anything as shown in Fig. 13. Therefore, we need to do some corrections for coupling multiple grains.

Step 3. Coupling of multiple level set functions. As in [15], we use the following correction for coupling multiple level set functions as

^i k;h = 1₂  ~i k;h min_j6=i ~k;hj  ; (6.11) where ^i k;h is the corrector of ~k;hi .

Step 4. Construct signed distance function by the reinitialization equation. In the step 2, the signed distance functions which are given initially do not satisfy the property of signed distance function any more. However, a number of simplications

can be made when i

k;hare signed distance functions. For this reason, a reinitialization

procedure is needed to reset the level set function to be a signed distance function. We discretize the reinitialization equation (5.13) by explicit scheme as

(^i k;h) ^k;hi
t = S(^k;hi )  1 (Dx^k;hi )2+ (Dy^k;hi )2 1 2 : (6.12)

For numerical purposes, a smooth sign function is used as in [15] S(^i k;h) = 8 > > > > < > > > > : 1 if ^i k;h > " ^i k;h " + 1  sin  ^i k;h " ! if j^i k;hj  " 1 if ^i k;h < " ; (6.13)

where " is a regularization parameter. Since the scheme is explicit, numerical stability requires a time step restriction as [15]

b 
t
x2 +
t
y2  < 1₂; (6.14)

where
x and
y are the grid sizes, b =  > 0 is a constant,  is the grain boundary mobility and  is the surface tension.

7 Numerical results

In this section, there are four numerical examples with Neumann boundary con-dition, and last one numerical example with periodic boundary condition. For the rst four examples, we consider the domain [ 0:5; 0:5]  [ 0:5; 0:5], and choose the mesh width
x =
y = 1=100, the constant b = 0:05, the parameter " =
x=2 and the time step
t = 0:001. In the last one example, we consider the domain [0; 3p3=5]  [0; 1:2], and choose the mesh width
x = (3p3=5)=100;
y = 1:2=100, the constant b = 0:025, the parameter " =
x=2 and the time step
t = 0:001.

7.1 Types of grain topological change

Let N denote the number of grains, V denote the number of vertices and E denote the number of grain boundaries. The grain structure in two-dimensional space obeys the following Euler equation [20]

N + V E = 1: (7.1)

The above geometry rule leads to a topological reconstruction of the entire grain structure. As in [21], there are two types of grain topological change. Fig. 14 shows the time evolution of grain topological change for type 1 by multiple level set method. At the initial time, we consider a triangle grain at the center of the domain. After the initial time, the shape of the central grain becomes like the three nodes cell. According to the motion by mean curvature, the shade of the central grain would shrink and even to a point. Fig. 15 shows the time evolution of grain topological

Figure 14: The time evolution of grain topological change for type 1 by multiple level set method.

Figure 15: The time evolution of grain topological change for type 2 by multiple level set method.

change for type 2 by multiple level set method. The grain topological change of type 2 involves the switching of the grain boundaries when two triple-junctions come very close to each other. Two grains lose an edge while two other grains gain an edge, thus maintaining the total number of edges and grains in the grain structure.

7.2 Von Neumann law

From the Von Neumann law, one can see that the grain with six edges will unchange. The grain with more than six edges will grow while the one with less than six edges will shrink. Now, we try to explain this motion from the geometric viewpoint. An equilateral hexagon has an interior angle with 120 degrees. Since the surface tensions are equal without surfactant, the angles at the triple-junction are same as the interior angle (120 degrees). Therefore, the grain with six edges will not change at all. A square has an interior angle with right angle. Since the surface tensions are equal without surfactant, the angles at the triple-junction must be expanded from right angle to 120 degrees. So each curve between two triple-junction points is convex. Due to the motion by mean curvature, the curves will at and the grain with four edges will shrink to a point as shown in Fig. 16. Fig. 17 shows the time evolution of eight nodes grain by multiple level set method. It is interesting to see that the central grain expands at beginning and shrinks later. Since an equilateral octagon has an interior angle more than 120 degrees, each curve between two

Figure 17: The time evolution of eight nodes grain by multiple level set method. junction points is concave. Due to the motion by mean curvature, the curves will at and the grain with four edges will expand. By using Neumann boundary conditions, the curves near the domain are perpendicular to the domain. When the wall number of the central grain becomes four, the central grain will shrink to a point.

7.3 Imposition of periodic boundary conditions

The above boundary condition is imposed by Neumann boundary condition. In this subsection, The level set equations are solved with periodic boundary condition. Firstly, we need that the initial signed distance functions are periodic on the domain boundaries. For clarity, we show that the contour of initial signed distance function for each grain as in Figs. 18 and 19. Excluding the boundary conditions, every numerical

Figure 18: The contour of initial signed distance functions with the periodic boundary conditions.

Figure 19: The contour of initial signed distance functions with the periodic boundary conditions.

Figure 20: The time evolution of grains with periodic boundary conditions by multiple level set method.

steps follow the previous section. The periodic boundary conditions are written by ~i m;h = ~1;hi ; 8 h = 1; 2; : : : m 1; ~i 0;h = ~m 1;hi ; 8 h = 1; 2; : : : m 1; ~i k;0 = ~k;m 1i ; 8 k = 0; 1; : : : m; ~i k;m = ~k;1i ; 8 k = 0; 1; : : : m: (7.2)

Fig. 20 shows the time evolution of grains with periodic boundary conditions by multiple level set method. The geometry includes two grains with seven edges, two gains with ve edges and eight grains with six edges. By von Neumann law, one can see that the gains with seven edges will expand and the gains with ve edges will shrink and disappear. In fact, the grain with more than six edges will grow and the one with less than six edges will shrink. After two gains with ve edges disappears, there are two grains with seven edges, one grain with eight edges, four grains with ve edges and three grains with six edges in the domain. While four grains with ve edges disappears, the domain includes two grains with seven edges, two grains with four edges and two grains with six edges. After two grains with four edges disappears, there are four grains with six edges in the domain only.

8 Conclusions

In this paper, we propose a nite dierence method to track curves motion whose normal velocity is determined by surface tension times the local mean curvature. We introduce the surfactant into the curves network and the surface tension varies fol-lowing surface diusion of the surfactant. The equations of motion are governed by a parabolic equation for the curve motion as well as a convection-diusion equation for the surfactant concentration along each curve. Our numerical method is based on direct discretization of the governing equations and the associated boundary con-ditions, which conserves the total surfactant mass in the curve network. Numerical examples are presented to illustrate how the inhomogeneous surface tension aects the motion of the curves and the evolution of the curve network. In the second part, we present a multiple level set method to track the grain topological change. The motion of grain boundary is considered by mean curvature. In this part, there is no surfactant on each grain boundary. The governing equations consist of a level set equation for the grain boundary motion and reinitialization equation for recon-structing signed distance function. Thus, one important step is coupling of multiple level set functions. Numerical examples shows the topological change evolution of the grain structure.

References

[1] L. Bronsard and B. T. R. Wetton, "A numerical method for tracking curve networks moving with curvature motion", J. Comput. Phys., 120 (1995) 66-87. [2] M.-C. Lai, Y.-H. Tseng and H. Huang, "An immersed boundary method for

interfacial ows with insoluble surfactant", J. Comput. Phys., 227 (2008) 7279-7293.

[3] B. Merriman, J.K. Bence and S. Osher, "Motion of multiple junctions - a level set approach", J. Comput. Phys., 112 (1994), 334-363.

[4] M. Sussman, P. Smereka and S. Osher, "A level set approach for computing solutions to incompressible two-phase ow", J. Comput. Phys., 114 (1994), 146-159.

[5] B. Engquist and S. Osher, "Stable and entropy satisfying approximations for transonic ow calculations", Math. Comp., 34 (1980), 45-75.

[6] P.G. DE Gennes, F. Brochard-Wyart and D. Quere, "Capillarity and Wetting Phenomena, Drop, Bubbles, Pearls, Waves", Springer (2004).

[7] M.P. Anderson, D.J. Srolovitz, G.S. Grest and P.S. Sahni, "Computer simulation of grain growth-I: kinetics", Acta Metall., 32 (1984), 783-791.

[8] P.S. Sahni, G.S. Grest and S.A. Safran, "Kinetics of the Q-state potts model in two dimensions", Phys. Rev. Lett., 50 (1983), 263-266.

[9] H. Huang, G.H. Gilmer and T.D. de la Rubia, "An atomistic simulator for thin lm deposition in three dimensions", J. Appl. Phys., 84 (1998), 3636-3649. [10] J.S. Chen, V. Kotta, H. Lu, D. Wang, D. Moldovan and D. Wolf, "A variational

formulation and a double-grid method for mesoscale modeling of stressed grain growth in polycrystalline materials", Comput. Methods. Appl. Mech. Eng., 193 (2004), 1277-1303.

[11] E. Fried and M.E. Gurtin, "A united treatment of evolving interfaces accounting for small deformations and atomic transport: grain-boundaries, phase transi-tions, epitaxy", TAM Reports, 1023 (2003).

[12] W.W. Mullins, "Tow-dimensional motion of idealized grain boundaries", J. Appl. Phys. 27, (1956) 900-904.

[13] W.W. Mullins, "in Metal Surfaces: Structure, Energetics, and Kinetics", W.D.Robertson and N.A.Gjostein, eds., American Society for Metals, Metals Park, Ohio, 1963, pp. 17-66.

[14] von Neumann, J. "in Metal Interfaces" (C.Herring, ed.) 108-110 (American So-ciety for Metals, Cleveland, 1951).

[15] X. Zhang, J.S. Chen and S. Osher, "A multiple level set method for modeling grain boundary evolution of polycrystalline materials", Interaction and Multi-scale Mechanics, 1(2), (2008) 191-209.

[16] Z. Pan, B. T. R. Wetton, "A numerical method for coupled surface and grain boundary motion", European J. Appl. Math. 19 (2008), no 3, 311-327.

[17] Y. Pawar, K.J. Stebe, "Marangoni eects on drop deformation in an extensional ow: the role of surfactant physical chemistry, I. Insoluble surfactants", Phys. Fluids 8 (1996) 1738.

[18] H. A. Stone, "A simple derivation of the time-dependent convective-diusion equation for surfactant transport along a deforming interface", Phys. Fluids A 2(1), (1990) 111-112.

[19] H. Wong, D. Rumschitzki and C. Maldarelli, "On the surfactant mass balance at a deforming uid interface", Phys. Fluids, 8(11) (1996) 3203-3204.

[20] C.S. Smith, "Grain shapes and other metallurgical applications of topology", Metal Interfaces, American Society for Metals, Cleveland, Ohio, (1952) 65-108. [21] J.E. Morral and M.F. Ashby, "Dislocated cellular structures", Acta Metall, 22