Regular article A level set based finite element algorithm ... - CiteSeerX

Comput Visual Sci 7: 97–110 (2004) Digital Object Identifier (DOI) 10.1007/s00791-004-0141-4

Computing and Visualization in Science

Regular article A level set based finite element algorithm for the simulation of dendritic growth Michael Fried∗ Institut für Angewandte Mathematik, Universität Freiburg, 79104 Freiburg, Germany (e-mail: [email protected], http://www.mathematik.uni-freiburg.de/IAM/homepages/micha) Received: 29 May 2001 / Accepted: 15 January 2003 Published online: 17 August 2004 –  Springer-Verlag 2004 Communicated by: G. Wittum

Abstract. Dendritic growth is a nonlinear process, which falls into the category of self-organizing pattern formation phenomena. It is of great practical importance, since it appears frequently and, in the case of alloys, affects the engineering properties of the resulting solid. We describe a new finite element algorithm for the two–dimensional Stefan problem, where the free boundary is represented as a level set. This allows to handle topological changes of the free boundary. The accuracy of the method is verified and several numerical simulations, including topological changes of the free boundary, are presented.

1 Introduction Dendritic growth is an important specimen of self–organizing pattern formation phenomena. Dendrites which arise during the solidification of a metallic melt affect engineering properties such as strength and durability of the later solid. Of less practical interest, but of fascinating beauty are the hexagonal structured and often intricated patterns of snowflakes, the dendrites of water. Here we consider the solidification of a pure, undercooled liquid. We base on a well–known generalization of the classical Stefan problem, cf. [2, 21, 23, 25, 38] . In order to fix notations, we recall the two–dimensional case of this model. Let Ω ⊂ R2 be a region filled with a pure substance, which may be in solid or liquid phase. The solid and liquid regions are then denoted by ωs and ωl , respectively. Together with the one–dimensional phase boundary γ they give a time-dependent partition p = (ωs , γ, ωl ) of the region Ω. This phase distribution p = p(t) and the dimensionless temperature θ = θ(x, t) satisfy the following evolution laws: ∗ Research supported by the Deutsche Forschungsgemeinschaft Schwerpunkt DANSE

θt − ∆θ = 0 in !! ""l ∂θ v=− on ∂n s β(n) v + a(n) k + θ = 0 on

ωs (t) ∪ ωl (t) ,

(1)

γ(t) ,

(2)

γ(t) ,

(3)

where n is the unit normal to γ into the liquid region ωl , v is the normal velocity of γ and k is the curvature of γ with a positive sign when ωs is convex. By [[·]]ls we denote the jump across the phase boundary γ with sign convention given by “liquid” – “solid”. The functions β and a denote material quantities: the kinetic coefficient β is defined on the unit sphere S, the interface potential a can be calculated from a smooth surface tension φ = φ(n) by the formula a(n) = D2 φ(n) n ⊥ · n ⊥ , where D2 φ(n) is the Hessian of the positively homogeneous extension of φ on R2 and ⊥ denotes the counterclockwise rotation by 90◦ . We assume that the temperature θ is scaled such that θ = 0 is the melting temperature of the regarded substance, that is the temperature of a planar phase boundary in equilibrium. In contrast to the classical Stefan problem, the above model includes a generalized Gibbs–Thomson law with undercooling, which allows for a shift in the temperature θ at the free boundary due to kinetic and curvature effects. Since the coefficients β and a depend on the normal n, the generalized Gibbs–Thomson equation also includes the modeling of directional anisotropies. In the sequel we suppose that the material quantities fulfill β∗ := inf β(n) > 0 n∈S

and a∗ := inf a(n) > 0 . n∈S

(4)

In the present article we consider the numerical solution of the two–dimensional modified Stefan problem (1), (2), (3) using a finite element method. The most difficulties for the numerical treatment of this Stefan problem arise from the generalized Gibbs–Thomson law (3). In view of this equation, the evolution of the free boundary γ(t) may be interpreted as

98

an anisotropic mean curvature flow with the temperature θ as a “driving force”. Therefore a common starting point to develop a numerical scheme for the dendritic growth consists of an algorithm for the mean curvature flow problem. Once such an algorithm is given, and if, on the other hand, we have a method to calculate the temperature θ provided the phase distribution p(t) is known, we may combine those methods to achieve a numerical scheme for the whole Stefan problem. We chose a numerical method for the mean curvature flow part of Problem (1)–(3) which uses a level set representation of the phase boundary γ(t) as it was first proposed as a computational procedure by Osher and Sethian for problems involving the curvature dependent evolution of an interface, cf. [27]. The basic idea behind such a level set method is to interpret the interface γ as the zero level of some continuous real valued function u. This enables us to handle topological changes of the phase boundary as they may occur for example if different parts of a dendrite are growing towards each other finally touching one another. In case of mean curvature flow like problems, the level set approach leads to a degenerate nonlinear partial differential equation for the function u, which has to be understood in a viscosity sense, cf. [3, 13–16, 20]. A broad variety of algorithms has been developed to numerically approximate the curvature dependent evolution of a hypersurface. A description of such methods could be found in [11] and the references therein. In Sect. 2 we outline a level set based finite element algorithm for the generalized mean curvature flow problem as it is given by the generalized Gibbs–Thomson law. Here we suppose for the moment the temperature θ to be known. Using a regularized formulation of the anisotropic and inhomogeneous level set problem, we combine time discretization by a semi implicit Euler method with the standard procedure to get a finite element method for the function u, and such for the generalized evolution of the free interface γ , which is defined by the zero level of u. Some numerical results are presented to illustrate the applicability of the algorithm for the anisotropic mean curvature flow of level sets. While an explicit calculation of the curvature of the level sets is not needed, in case of several applications it is desirable to know the curvature of the moving interface γ . Particularly with regard to our algorithm for the modified Stefan problem, where we have to know the curvature explicitly, we end Sect. 2 by briefly describing a numerical method to approximate the mean curvature of a given finite element function’s level sets. Vice versa we assume in Sect. 3 the phase boundary γ(t), its normal velocity v and its curvature k to be known. Then a numerical method for approximating the second unknown, the temperature, is obtained by a similar finite element method as it was introduced by Schmidt in [30] and [31]. The difference lies in the fact that Schmidt used a parametric ansatz for the representation of the phase boundary, while we represent γ using a level set method. To obtain a numerical method for the full modified Stefan problem we combine in Sect. 4 the finite element methods for the free boundary and for the temperature. We end by presenting various numerical results obtained by our method. They include a verification of its accuracy with the help of an example with explicitly known exact solution and an example with topological changes of the phase boundary. The described algorithms were implemented with

M. Fried

the help of ALBERTA, an adaptive multi–level finite element toolbox which uses bisectioning refinement and error control based on residual techniques, and provides also a wide range of linear system solvers, cf. [32]. 2 A level set formulation for the mean curvature flow In this section we suppose the temperature θ to be given. Then the generalized Gibbs–Thomson law (3), leads to the following generalized mean curvature flow for the evolution of the phase boundary γ(t): Problem 1. Find a family of interfaces {γ(t)}t∈[0,T ) such that β(n)v + a(n)k = −θ

on

γ(t), t ∈ [0, T ) ,

(5)

where γ(0) = γ0 . The purpose of this section is to develop a finite element algorithm for the anisotropic and inhomogeneous mean curvature flow given by the above problem, which can be extended to an algorithm for the Stefan Problem (1)–(3). As mentioned in the introduction, we consider a level set approach to the mean curvature flow Problem 1, assuming the free boundary γ to be represented by the zero level of some continuous real valued function u. In addition we assume that all level sets of u evolve according to equation (5). A formal calculation then leads to the following non–linear degenerate and singular parabolic partial differential equation for the function u, compare [13]: # $ # $ ∇u ut ∇u ∇u β −a ∇· =θ. (6) |∇u| |∇u| |∇u| |∇u| Conversely, if u is a viscosity solution of (6) in R2 , then each level set of u evolves in a weak sense as it is required by (5). Existence and uniqueness of a viscosity solution u of (5) with an initial condition u(·, 0) = u 0 ∈ C(R2 ) in arbitrary spatial dimensions ≥ 2 was shown by Evans and Spruck, who investigated the isotropic and homogeneous problem in a series of publications ([13–16]). At the same time, similar results were obtained by Chen/Giga/Goto and Giga/Goto/Ishii, even for anisotropic and inhomogeneous situations, cf. [3, 20]. In view of Problem 1, the above authors showed that not only the viscosity solution u is uniquely determined by (6) with the above initial condition, but also the generalized evolution of γ(t), defined by {x ∈ R2 |u(x, t) = 0}, does only depend on initial hypersurface γ0 = γ(0), and therefore is independent of the special choice of the initial value u 0 . These results allow to reformulate Problem 1 in a level set context: Problem 2. Given some function u 0 with γ(0) = {x ∈ R2 | u 0 (x) = 0} , find the continuous function u : R2 × [0, T ) −→ R such that # $ # $ ∇u ut ∇u ∇u β −a ∇· =θ, |∇u| |∇u| |∇u| |∇u| with initial condition u(·, 0) = u 0

∈ R2 .

A level set based finite element algorithm for the simulation of dendritic growth

There are various methods to numerically solve Problem 2, see for example [5, 27, 33, 39]. One of the main difficulties in discretizing this problem arises if the gradient ∇u vanishes. To avoid such difficulties, we turn to the following regularization of Problem 2: Problem 3. For fixed ε > 0, find u ε ∈ C 0 (R2 × [0, T )) such that βε (∇u ε )

u εt ∇u ε − aε (∇u ε ) ∇ · =θ Q(u ε ) Q(u ε )

(7)

in R2 ,

where Q(u ε ) is given by % Q(u ε ) = ε2 + |∇u ε |2

(8)

(9)

The functions βε and aε denote smooth extensions βε , aε ∈ C ∞ (R2 ) of the given anisotropies β and a, respectively:  # $ p   , | p| ≥ ε β   # | p| # $$ # $ 2 p p 1− 2 ε 2 βε ( p) := ε −| p| + β β − β e , 0 < | p| < ε  0  | p| | p|    β0 , | p| = 0 , with β0 = 12 (minS β + maxS β), and the analogous definition for aε . Suppose β, a ∈ C ∞ (S) fulfilling (4) and u 0 ∈ C ∞ (R2 ) is constant outside of some large ball, then for every 0 < ε < 1 there exists a unique bounded solution u ε of Problem 3 with sup -u ε , ∇u ε , u ε,t -L ∞ (R2 ×[0,T )) ≤ C,

(10)

0 0 in ωl (0), u 0 < 0 in ωs (0), where ωs (0) denotes the part of the domain Ω which lies inside of the initial interface γ0 , and ωl (0) := Ω \ (ωs (0) ∪ γ0 ). Remark 1. Problem 3 is closely related to the mean curvature flow of a graph, compare i.e. [24, 26, 37]. 2.1 A finite element scheme for the mean curvature flow In order to develop a finite element algorithm, we turn now to the following weak formulation of the regularized problem: find u ε (·, t) ∈ H 1 (Ω) such that for all φ ∈ H01 (Ω) * * * βε (∇u ε ) u ε,t φ ∇u ε ∇φ θφ + = , (13) aε (∇u ε ) Q(u ε ) Q(u ε ) aε (∇u ε )

Ω

Ω

Ω

with initial conditions (8) and Dirichlet boundary data. Let T be a conforming non-degenerate simplicial triangulation of the polygonal domain Ω. We denote by Pk the space of polynomials of degree ≤ k and define the finite element spaces V by + , . V := Φ ∈ C 0 Ω | ∀S ∈ T , Φ| S ∈ P2 (S) , ◦

and V := V ∩ H01(Ω). For simplicity we chose a fixed timestep size/ τ0> 0. Using the notations t m = mτ for m = 0, 1, . . . , M = Tτ and f m = f(t m ) for a given time dependent function f , a semi– implicit scheme for the anisotropic mean curvature flow of level sets is given by: Algorithm 1. For given U0 (·, t) ∈ V, t ∈ [0, T ), find U m ∈ V (m = 1, . . . , M) such that * m−1 , m * U − U m−1 1 βε ∇U m ∇Φ Φ + τ aεm−1 Q m−1 Q m−1 Ω

Ω

=

*

Ω

m

θ Φ aεm−1

◦

∀Φ ∈ V ,

(14)

100

M. Fried

with initial and boundary conditions ◦

U 0 = U0 , U m (·) − U0 (·, t m ) ∈ V , where U0 (·, t) ∈ V denotes a finite element representation of u 0 (·, t). In contrast to the highly nonlinear implicit scheme obtained by using Q m instead of Q m−1 in (14), in case of the above semi–implicit method we only have to solve a linear problem per timestep, while on the other hand a restriction on the smallness of the timestep size with respect to the gridsize as it is necessary for the explicit scheme (∇U m−1 instead of ∇U m ) is not needed. Remark 2. In case of the isotropic and homogeneous situation, Deckelnick and Dziuk proved in [7] and [8] convergence results for algorithm (14) using piecewise linear finite elements. An investigation of the corresponding implicit, semi–implicit and explicit schemes for the mean curvature flow of a graph could be found in [10]. 2.2 Numerical experiments The following numerical experiments may illustrate the capabilities of our finite element algorithm for the isotropic as well as for the anisotropic mean curvature flow of level sets. The first example deals with an isotropic and homogeneous situation, where the exact solution is known. This enables us to verify the algorithm by calculating the experimental order of convergence. By the second experiment, we again investigate an isotropic and homogeneous case, to show the applicability of our scheme in situations where γ changes its topological type. The level set method not only covers topological changes by breaking or merging of the moving interface but also situations where the free boundary develops a nonempty interior, we show this by computing an inhomogeneous example introduced by Paolini and Verdi, cf. [29]. Using a threefold anisotropy, we approximate in the forth experiment the anisotropic mean curvature flow of a circle toward the Wulff shape of the given anisotropy. As an outlook we sketch in the last example the applicability of our algorithm in case of three dimensional situations like the shrinking of a torus. The shrinking circle. A known solution of the homogeneous and isotropic mean curvature flow is the shrinking 1 sphere γr0 (t) := {x ∈ R2 ||x| = r02 − 2t} with initial radius r0 > 0. For r0 = 1, the domain Ω = [−3.5, 3.5]2 and the initial function  2π 3 2π 3  r 2 + cos for r ≤ 3 , − cos 9 9 u 0 (r) = 2π 3   1 + cos for r > 3 , 9

we have the radial symmetric solution  2π , 2π 3 -3 ,  r 2 + 2t + cos for r 2 + 2t ≤ 9 , − cos 9 9 u(r, t) = 2π 3   1 + cos else , 9 where r = |x|. Here we choose Neumann boundary conditions and calculate the piecewise quadratic finite element

approximation U using several uniform triangulations T j . Starting from a given triangulation T 0 , the triangulation T j is obtained from T j−1 by twice bisecting each triangle in T j−1 , hence the gridsize h j is given by h j = 12 h j−1 . We performed several numerical tests, choosing the timestep size to be τ = 0.1 ∗ h, τ = h 2 and τ = h 3 , comparing different ε–h relations. Tables 1–3 show the rates of convergence obtained using different relations between the mesh size h j and the regularization parameter ε. The experimental order of convergence is defined by 2 3 Err ln Err j j−1 EOC j := , (15) ln(2) with error Err j := max -u(·, t m ) − U m -L 2 (Ω) . 1≤m≤M

For ε = h, we observe first order or lower rates of convergence with respect to h, while ε = h 2 and ε = h 3 lead to significantly higher rates of convergence in the range of 1.5 up to 2.5, at least using τ = h 2 or τ = h 3 . Even if the highest experimental order of convergence was reached using the combination ε = h 3 and τ = h 3 , for practical computations with fixed uniform grids we decided to choose τ = h 2 and ε = h 2 in order to achieve an efficient algorithm: choosing τ = h 3 results in a very small timestep size Table 1. EOC for ε = h, ε = h 2 and ε = h 3 . Time step size τ = 0.1h h

Err ε=h

EOC ε=h

Err ε = h2

EOC ε = h2

Err ε = h3

EOC ε = h3

1.75 0.88 0.44 0.22 0.11 0.05

1.357 0.890 0.631 0.383 0.207 0.112

– 0.61 0.50 0.72 0.88 0.89

1.428 0.840 0.358 0.129 0.052 0.026

– 0.76 1.23 1.47 1.32 1.00

1.455 0.789 0.212 0.076 0.040 0.024

– 0.88 1.90 1.49 0.93 0.76

Table 2. EOC for ε = h, ε = h 2 and ε = h 3 . Time step size τ = h 2 h

Err ε=h

EOC ε=h

Err ε = h2

EOC ε = h2

Err ε = h3

EOC ε = h3

1.75 0.88 0.44 0.22 0.11 0.05

3.945 1.674 0.841 0.425 0.208 0.109

– 1.24 0.99 0.99 1.03 0.94

3.987 1.628 0.578 0.179 0.054 0.018

– 1.29 1.49 1.69 1.73 1.62

4.005 1.581 0.454 0.133 0.043 0.015

– 1.34 1.80 1.77 1.68 1.59

Table 3. EOC for ε = h, ε = h 2 , and ε = h 3 . Time step size τ = h 3 h

Err ε=h

EOC ε=h

Err ε = h2

EOC ε = h2

Err ε = h3

EOC ε = h3

1.75 0.88 0.44 0.22 0.11 0.05

4.079 1.536 0.600 0.364 0.199 0.105

– 1.41 1.36 0.72 0.87 0.92

4.136 1.489 0.363 0.108 0.028 0.007

– 1.47 2.03 1.76 1.96 1.99

4.160 1.441 0.241 0.048 0.009 0.001

– 1.53 2.58 2.33 2.48 2.54

A level set based finite element algorithm for the simulation of dendritic growth

hence in a long overall computation time (i.e. the computation time for the single run of the last row of Table 3 was more than 36 h, while the corresponding computation for Table 2 took just about 4h). Similar, choosing a small ε usually leads to a higher number of cg iterations to solve the linear system, i.e. ε = h 2 leads to an average of 34 iterations compared to an average of 39 iterations in case of ε = h 3 , hence again to a longer computation time. Topological changes. One advantage of the level set approach is the ability to handle topological changes during the evolution of the free boundary. To illustrate this feature of our algorithm, we approximate in the following experiment the mean curvature flow of a lemniscate. We consider γ0 to be given by the zero level of the function -2 , 2 , x 1 + (a x 2 )2 − 2 x 12 − (a x 2 )2 u 0 (x) = , 30 000 with x = (x 1 , x 2 ) ∈ R2 and a scaling parameter a ∈ R+ . The zero level of this function looks like an “∞”, compare Fig. 1. Depending on the initial choice of the parameter a ∈ R+ , we found two qualitatively different types of topological changes: If a > 1, the initial lemniscate γ0 is thin and long, thus the outer angles α1 are bigger than the inner ones, cf. Fig. 1. Given such initial data, the evolving curve breaks into two tear shaped parts, which independently of each other shrink to two points before vanishing. Figure 2 shows the evolution in the case a = 1.5. If, on the other hand, we choose the parameter a to be a < 1, then the angles α and α1 of the initial lemniscate are α1 < α, and the interface breaks at the self intersection such, that it stays one connected curve which shrinks down to one point at the origin, cf. Fig. 3. Fattening. Evans and Spruck conjectured that instead of breaking, an evolving interface may develop a nonempty interior, if at some time the interface is not the boundary of a smooth open set, cf. [13, 15]. This phenomenon is called fattening. Simply speaking, fattening is the “getting flat” of the function u around its zero level. In our numerical experiments concerning the evolution of a given interface, we always checked for the appearance of such an effect in a rather rough way by testing if the finite element function U vanishes on some simplices of the triangulation. Since this did not promise very strict results, in the next example we used a more accurate method to numerically check for fattening. It was invented by Paolini and Verdi, cf. [29], and is based on the following idea: If the zero level of the function U develops a nonempty interior, one will find a jump in the volume At (λ) := vol({x ∈ Ω|U(x, t) ≤ λ}) at λ = 0. Paolini and Verdi

Fig. 1. Lemniscate

101

Fig. 2. Shrinking lemniscate with parameter a = 1.5

Fig. 3. Shrinking lemniscate with parameter a = 0.5

calculated this quantity to show the numerical evidence of fattening using a phase field approach to approximate the inhomogeneous mean curvature flow of two circles. At certain critical time during the evolution, the two circles are touching each other. We investigated this example by the help of our level set based algorithm. In contrast to the previous examples we are now in a situation with non vanishing right hand side θ. Paolini and Verdi considered the mean curvature flow of two circles with centers at (x m , 0) and (−x m , 0), and initial radius r0 . Choosing the function θ = θ(t) := 2 − t, the evolution of each of those circles isolated from the other one is given by the ordinary differential equation r 1 (t) +

1 = θ(t) , r(t)

r(0) = r0

for the radius r(t). For r0 smaller than a critical radius rcrit , also the evolution of both circles is ruled by the above or1 dinary differential equation. Until the time t = 2 − r(t) the radius r(t) is growing, while for later times it shrinks until the circles vanish. If the initial radius r0 > rcrit , then there is a time t∗ at which r(t∗ ) = x m holds. In this case, the evolution of one isolated circle differs from that of two circles, which will touch each other and merge to one single connected curve. At the time t∗ , we are in a similar situation as we were at the initial situation of the previous example. Due to the results of Paolini and Verdi, we expect that the interface will fatten up, developing a nonempty interior. Calculating the volumina At m (λ) for λ ∈ [−0.01, 0.01] at a time step t m with |t∗ − t m | ≤ τ, we indeed found the numerical evidence for fattening of a level set which evolves by its mean curvature, cf. Fig. 4. This result corresponds to the results obtained by Paolini and Verdi, cf. [29]. As fattening of the zero level of the computed u means that the gradient ∇u vanishes (or at

102

M. Fried

20 19 18 17 16 15 14

Fig. 5. c = 0.8: Wulff shape of the anisotropy a

13 12 -0.01 -0.008-0.006-0.004-0.002 0 0.002 0.004 0.006 0.008 0.01 Fig. 4. Fattening. Volumina A0.45 (λ) for −0.01 ≤ λ ≤ 0.01

least becomes very small) around the zero level set of u, the influence of the regularization parameter ε in the regularized % modulus of ∇u, Q(u) = |∇u|2 + ε2, becomes more important. Therefor, we did similar computations using different values of the regularization parameter ε = 10−1 , 10−3 , 10−7 . In all three cases, we found the same behavior. Nevertheless, it should be noted, that this result does not necessarily imply a “real” fattening of the function u which might be very flat around a still one dimensional zero level line. Also this does not imply a fattening of the continuous viscosity solution to that problem. The described result should be taken as a hint that during the computation of the regularized mean curvature flow there might occur a fatten up of the discrete free boundary. Anisotropic mean curvature flow. The following example illustrates the application of our algorithm in case of the anisotropic mean curvature flow. To keep it simple, we suppose β to be constant β = 1 and a to be given as a(n(α)) = f(α) + f 11 (α) = 1 + c · cos(3α) ,

where α denotes the angle between the normal n and the x 1 – axis, c is a given parameter c > 0 and f(α) = 1 − 8c cos(3α). A method to visualize the amount of anisotropy imposed by the function a is given by depicting the the Wulff shape + . W a := ∂ p ∈ R2 |2 p, n(α)3 ≤ f(α) ∀α ∈ R .

In some sense, the Wulff shape W a can be viewed as the unit sphere of the anisotropic metric given by the anisotropy a. In the isotropic situation c = 0, the Wulff shape W a fulfills W a = ∂B1 (0) . Due to the assumption (4) the anisotropy a is convex in the sense of Gurtin, cf. [22], and in case of a convex anisotropy it is known that a convex initial curve γ0 during the evolution approximates the Wulff shape W a , cf. [19, 28]. For the following anisotropic experiment, we choose c = 0.8. Figure 5 shows the Wulff shape with threefold symmetry belonging to that choice of the parameter c. We start with the unit circle and calculate the anisotropic evolution on the domain Ω = [−4, 4]2 choosing Neumann boundary conditions at ∂Ω. The shrinking circle converges to the Wulff shape rather quickly, compare Fig. 6.

Fig. 6. c = 0.8: Anisotropic shrinking circle

2.3 Computing the curvature Using the level set method to approximate the mean curvature flow of a given interface, an explicit calculation of the interface’s mean curvature is not needed. However, for certain applications – as for instance for the modified Stefan problem – one would like to know the curvature of γ explicitly, as it appears in the computation of some other quantities, compare Sect. 4. Unfortunately, in our situation there is no straightforward definition of the curvature of γ(t m ), as the free boundary is given by a level set of a piecewise polynomial function U m ∈ V . One possible letout is the definition of a discrete curvature K ∈ V of the level sets of a given function U ∈ V by * * ◦ ∇U m ∇Φ KΦ = − ∀ Φ ∈ V, (16) m Q(U )

Ω

Ω

with appropriate boundary values K = K 0 on ∂Ω. This definition is motivated by the way of regularizing the mean curvature equation (6) and the well known formula for the curvature of the level sets of a smooth function u with non ∇u vanishing gradient k = ∇ · |∇u| . In case of the below presented numerical experiments, we chose K 0 to be a finite element approximation of k0 (x) = |x|−1 , which represents the curvature of the level

A level set based finite element algorithm for the simulation of dendritic growth

sets of a rotationally symmetric function. In principle, U m is not a rotationally symmetric function, but a discrete distance function from γ(t m ) (compare Sect. 4.2 below), and so the level sets of U m near the boundary ∂Ω are approximately circular, at least as long as the interface γ(t m ) is far away from ∂Ω. Numerical examples. The following numerical experiments illustrate numerical convergence as well as possible problems of the above finite element approximation of the curvature. Firstly we consider the function u(x) = − cos( π9 |x|2 ) + π9 . As the function u is rotational sym1 . metric, the curvature of its level sets is given by k1 (x) = |x| We calculate the finite element approximation K 1 of the interpolation U ∈ V of u on the domain Ω = B2 (0) and compute the experimental order of convergence EOC using linear, quadratic, and cubic finite elements. Again, we define the EOC by (15) using L 2 -errors. Since the curvature k1 has a singularity at the origin, we compute the L 2 -error on the domain D = Ω \ B0.5 (0). For all computations we fixed the parameter ε to be ε = 0.0001. Table 4 shows the obtained results. Obviously using linear finite elements the method did not converge, we need to choose at least piecewise quadratic elements for k as well as for the discrete level set function U. Further problems may arise if the gradient of the function vanishes on a whole level line γc ∈ Ω. We choose the function 4 −dist(x, ∂B1 (0) ∪ ∂B2(0)) for 1 ≤ |x| ≤ 2 d(x) := dist(x, ∂B1 (0) ∪ ∂B2(0)) else , on the domain Ω = B2 (0). Again the function is a rotationally symmetric, but this time the gradient ∇d vanishes on the level line d = 12 in between the both parts of the zero level set. Evaluating the L 2 -errors on the domain D as before, even for cubic elements we obtain no convergence of the finite element approximation K 2 to the curvature k2 of the level lines of d. But the error could be localized around the level line with vanishing gradient, by omitting the ring R := B1.75(0) \ B1.25 (0) and calculating the error on the domain D1 := D \ R we found the results shown in Table 5. Table 4. L 2 -errors and EOCs for the approximation of k1 hj 1 0.5 0.25 0.125 0.0625 0.03125

Err j linear

EOC j linear

Err j quadratic.

EOCj quadratic.

7.6639 9.6335 12.0309 12.4556 12.6097 12.6950

– −0.33 −0.32 −0.05 −0.02 −0.01

2.3349 0.9301 0.4652 0.2428 0.1228 0.0615

– 1.33 0.99 0.94 0.98 0.99

Table 5. L 2 -errors and EOCs for the approximation of k2 hj 1.0 0.5 0.25 0.125 0.0625

Err j in D

EOC j in D

Err j in D1

23.0293 28.4363 45.4699 54.8327 78.9481

– −0.30 −0.68 −0.27 −0.53

2.4 2.5759 0.5952 0.2482 0.1736

EOCj in D1 – −0.10 2.11 1.26 0.52

103

3 Finite element method for the heat equation The previous section dealt with the anisotropic mean curvature flow of level sets, assuming the right hand side θ to be known. Vice versa, in this section we suppose that the free boundary γ(t), its normal velocity v and curvature k are given. Then, together with boundary and initial conditions, the equations (1)–(3) determine the temperature θ. Problem 4. Find the temperature θ such that θt − ∆θ = 0 v=−

!!

∂θ ∂n

in [Ω \ γ(t)] × (0, T ) ,

""l

on γ(t) ,

(18)

s

β(n) v + a(n) k + θ = 0

on γ(t) ,

(19)

/ 0 on Ω × {0} ∪ [∂Ω × (0, T )]

θ = θ0

(17)

(20)

In order to derive a weak formulation of Problem 4, we pursue the method proposed by Schmidt, cf. [31]. After integration by parts and utilizing the equations (18), (19) we get for all φ ∈ H01(Ω) and all t ∈ (0, T ) * * * * 1 a θt φ + ∇θ ∇φ + θφ=− kφ (21) β β

Ω

Ω

γ(t)

γ(t)

The spatial discretization is done using the finite element ◦ spaces V and V as defined in Sect. 2.1 and an fully implicit scheme in time: Algorithm 2. (Schmidt, 1996) Given smooth interfaces γ(t m ) with curvature k and normal n and boundary values Θ0 (·, t) ∈ V, t ∈ [0, T ), find Θ m ∈ V (m = 1, . . . , M) such that * * * 1 , m 1 m Θ − Θ m−1 Φ + ∇Θ m ∇Φ + Θ Φ τ β Ω Ω γ(t m ) * ◦ a(n) =− kΦ ∀Φ ∈ V , (22) β(n) γ(t m )

with initial and boundary conditions Θ 0 = Θ0 ,

◦

Θ m (·) − Θ0 (·, t m ) ∈ V ,

where Θ0 (·, t) ∈ V denotes a finite element interpolation of θ0 (·, t). A combination of this numerical scheme for the heat equation and the scheme 1 for the mean curvature flow leads to an finite element algorithm for the modified Stefan problem. 4 Finite element method for the Stefan problem We now turn back again to the full Stefan problem (1)–(3). In contrast to the assumptions of the previous sections, neither temperature nor phase boundary and it’s curvature are known.

104

M. Fried

In order to compute this quantities, we combine the described methods for the heat equation, the interface motion and the curvature. 4.1 Combining the algorithms Suppose for the finite element approximations Θ m−1 and U m−1 for the temperature θ(t m−1 ) and the level set function u(t m−1 ), respectively, to be known. Then the new values at the time t m are computed as follows: First, we solve equation (14) for U m using Θ m−1 in lieu of θ(t m ). The approximation γ m of the new phase boundary then is given by the zero level set of U m . Second, following the method described in Sect. 2.3, we calculate the discrete curvature K ∈ V . The time step t m ends by computing Θ m according to * * * 1 , m 1 m Θ − Θ m−1 Φ + ∇Θ m ∇Φ + Θ Φ τ β Ω

γm

Ω

=−

*

γm

a KΦ β

◦

∀Φ ∈ V ,

(23)

which is a modified version of (22) using γ m and K in place of γ(t m ) and k, respectively. From Sect. 2.3 we know, that replacing the curvature k by the discrete curvature K ∈ V forces us to choose at least quadratic finite elements. Summarizing, this gives the following level set based finite element discretization of the Stefan problem (1)–(3): Algorithm 3. For given U0 (·, t) and Θ0 (·, t), find U m and Θ m ∈ V (m = 1, . . . , M) such that , * * βε ∇U m−1 U m − U m−1 1 ∇U m ∇Φ , , - Φ+ , τ aε ∇U m−1 Q U m−1 Q U m−1 Ω

=

Ω

*

Θ m−1 , -Φ aε ∇U m−1

Ω

◦

∀Φ ∈ V ,

(24)

and * * * 1 , m 1 m m−1 m Θ −Θ Φ + ∇Θ ∇Φ + Θ Φ τ βε Ω

γm

Ω

=−

*

γm

aε KΦ βε

◦

∀Φ ∈ V ,

(25)

with initial and boundary conditions ◦

U 0 (·) = U0 (·, 0) ,

U m (·) − U0 (·, t m ) ∈ V ,

Θ 0 (· = Θ0 (·, 0) ,

Θ m (·) − Θ0 (·, t m ) ∈ V ,

where

◦

+ . γ m := x ∈ Ω|U m (x) = 0 , and K ∈ V is given by * * ∇U m ∇Φ KΦ = Q(U m )

Ω

Ω

with boundary condition

◦

∀Φ ∈ V,

◦

K − K0 ∈ V . 4.2 Numerical algorithms We implemented the finite Element Algorithm 3 utilizing the toolbox ALBERTA, which provides adaptive multi–level finite element methods using bisectioning refinement and error control by residual techniques, compare [32]. In order to apply the discrete curvature defined in Sect. 2.3, we chose quadratic finite elements for all in the following described numerical experiments. 4.2.1 Approximating the interface. Possibly the most interesting quantity of the Stefan Problem, the phase boundary γ m , is implicitly given as the zero level set of the piecewise quadratic function U m , hence it is a piecewise conic section. Though the analytical determination of the levelset γ m is straightforward, the numerical determination is costly. In practice, we therefore use a polygonal approximation Γ m of γ m . Our first ansatz for such an approximation is to compute the zero levelset of a piecewise linear interpolation U˜ m m of U m on the regular refined grid T˜ . The zero levelset of m U˜ can easily be computed by calculating its zeros on the m edges of the triangles of the finer triangulation T˜ and conm necting them. Provided that U˜ is not vanishing on a whole triangle, this leads to a polygonal approximation Γ˜ m of γ m . Note that neither the polygon Γ˜ m has necessarily any common point with the zero levelset of U m , nor the topological type of γ m has to be preserved by Γ˜ m . In practice, this ‘linear’ method turned out to be not accurate enough to allow a stable evolution of the discrete phase boundary. On the left hand side of Fig. 7 a typical situation of an isotropic calculation is shown. A few timesteps after starting the calculation with a non-convex initial interface, instabilities occur and parts of the solid phase detach. To achieve a stable evolution of the interface using the above ‘linear’ method, we have to reduce the timestep size from τ = 0.01 to τ ≤ 0.00001, thus having thousandfold the numerical effort as before, which prohibits the use of this method for our purpose. One idea to improve the approximation of the interface is as follows: Instead of determing the zeros of the piecewise m linear interpolation U˜ m on the edges of the refined grid T˜ , m we now search for the zeros of the original function U on the edges of the triangulation T m . Since U m is a piecewise quadratic function, we may find up to two separated zeros on a single edge and hence connecting the zeros is not as simple as in case of the linear interpolation U˜ m . To connect them to achieve a polygonal approximation Γ m while preserving the topological type of the level set γ m , we apply an element wise working algorithm, which combines information given by the special distribution of the zeros on the triangles’ edges and the direction of gradient ∇U m at the zeros as well as it calculates additional zeros on certain straight lines inside of the triangles. Thereby we neglect such parts of the interface γ m , which are ellipses laying completely inside of a single simplex. A detailed description of this algorithm goes beyond the scope of this article, it can be found in [17]. Numerical experiments indicate that using this ‘quadratic’ method noticeable betters the stability of the interface’s evolution, compare the right picture of Fig. 7. Here we use the same parameters as we did before including the timestep size τ = 0.01. The discrete

A level set based finite element algorithm for the simulation of dendritic growth

105

where the sum is over all edges e of the triangle T , which belong to the interior of Ω, and the global error indicator ; κ ηmax mark T for refinement sum := sum + η2T end if end if end for end while A similar method leads to the subset B of triangles, which could be coarsened, compare [9]. Finally, the modification of the mesh is done by a bisectioning algorithm, which is taken from [1]. The initial triangulation is generated by an application of this adaptive procedure to the stationary problem * * * * 1 0 a ∇Θ 0 ∇Φ + Θ Φ = − Θt0 Φ − K Φ , (26) β β

Ω

Γ(0)

Ω

Γ(0)

◦

for Φ ∈ V with given initial data Γ(0) and Θt0 and boundary value Θ0 .This method was proposed by Schmidt, [30]. Beside the generation of an initial adapted grid, defining Θ 0 as the solution of (26) ensures that the initial temperature is compatible with the discrete interface Γ(0) and its curvature K , compare [30]. The finite element algorithms for the approximation of the regularized level set function u ε are analogous to the methods described for the heat equation. Namely, in both problems we apply the same finite element spaces, and consequently the same triangulation. Hence the adaptive process is completely done using the error indicator for the heat equation presented above. 4.2.3 Redistancing. Remembering the definition of the discrete curvature K , we see that this quantity depends on the

106

M. Fried

gradient ∇U m of the level set function U m . In practice, very steep or flat gradients ∇U m would lead to inaccurate approxi∇U m mations of the curvature and the normal n = -∇U m - , as well as this could cause stability problems like tip splitting. Therefore it is desirable to prevent having steep or flat gradients develop in U m , whenever this would be possible. For instance, if U m is a distance function to γ(t m ), we will avoid such ugly gradients. However even if we choose u 0 to be the signed distance function 4 dist(x, γ0) if x ∈ ωs (0) dγ0 (x) := x ∈ ωl (0) −dist(x, γ0) if from γ0 , the level set function will cease to be a distance function after one timestep. This may be fixed by reinitializing the function U m to be a distance function from the discrete phase boundary γ(t m ) at each timestep. There are several possibilities of redistancing the level set function, see for example Chopp [4], Strain [35] or Sussman, Sekerka and Osher [36]. In view of the fact that we already know a polygonal approximation Γm of the interface γ(t m ), we use a direct method, calculating a Lagrange interpolant of the distance function from the polygonal curve Γm . 4.3 Numerical results We implemented a version of the above described algorithm using the finite element tool-box ALBERTA, which was developed by Schmidt and Siebert, [32]. Here, we present numerical results obtained by this implementation. 4.3.1 Convergence of the numerical method. As a first test we compare the solution of our finite element method with the known solution in case of an isotropic test problem, compare [30]. This enables us to check the accuracy of our numerical method. Using the notations 1 2β ˙ f(t):= W(t) = R(t) := R02 + 2t , R0 ∈ R+ , R(t)3 * s − z2 √ 2β e 2 + W(t) := − , β∈R , T(s):=− e dz . R(t) z 1

the problem reads as follows  θt − ∆θ = f for t ≥ 0, x ∈ R2 \ Γ(t) ,     ! ∂θ " + VΓ = 0 on Γ(t) , ∂νΓ Γ(t)      θ + βCΓ + βVΓ = 0 on Γ(t) .

The difference between the modified Stefan problem (1)–(3) and the above problem is the right hand side f in the heat equation, which is chosen in such a way, that an explicit solution is known. This solution is given by  for |x| ≤ R(t)  W(t) # $ θ(x, t) := |x|  else W(t) + T R(t)

Table 6. Exact solution: growing circle. evolution of the error while uniformly refining the grid h ∆t E ∞,2 (Θh )

1.0 0.25 0.8319

0.5 0.0625 0.3546

0.25 0.0156 0.0657

0.125 0.0039 0.0256

0.0625 0.00098 0.0133

Fig. 8. Evolution of the radius using different uniformly refined grids

and the free boundary Γ(t) = {x ∈ R2 | -x- = R(t) }, where 1 1 v(t) = Rt (t) = R(t) and k(t) = R(t) are the normal velocity and curvature of the free boundary. Using the parameter β = 0.1 and the regularization parameter ε = 10−6, we investigate the behavior of error between the computed temperature Θ and the exact temperature θ, solving the same problem with different globally refined spatial meshes and time step sizes, c. f. [18]. The error was measured by = = E ∞,2 (Θh ) := max = Θ(·, ti ) − Θhi = L 2 (Ω) . i=0,...,n

Table 6 shows the results of these calculations. The evolution in time of the radii of Γh for this experiments is depicted in Fig. 8, here the solid line shows the exact evolution of the radius R(t).

4.3.2 Numerical experiments. In this section, we show the results of different numerical experiments in the case of the Stefan problem (1)–(3), where an explicit solution is not known. For all of the following experiments, we choose the boundary values Θ0 (·, t) = T0 to be constant on ∂Ω for all t ∈ [0, T ], and the functions a and β of the form β(n) = a(n) = δ(1 + A cos(kα + α0 )) , where α = αn denotes the angle between n and e1 -axis. The initial temperature Θ 0 is computed from (26) with Θt0 = 0. The first numerical experiment is a further isotropic calculation choosing the parameter A = 0. As opposed to the above test problem, we investigate the evolution of an initially non– convex phase boundary, namely > # $ ? cos(α) Γ(0) = a(n) | α ∈ [0, 2π] , (27) sin(α)

A level set based finite element algorithm for the simulation of dendritic growth

where a(n) = 0.15 + 0.4 cos(4α + π4 ). The investigation of such a curve was proposed by Sethian and Strain in [34]. Here, we are interested in the influence of the parameter δ. In some sense this parameter plays the role of the anisotropic functions β(n) and a(n), and in a realistic setting, it could be a very small parameter. We compare the numerical results for different values of the parameter δ. Figure 9 shows the evolution for the parameters δ = 0.006 and δ = 0.0009 respectively. For larger values of δ as in the first case (Fig. 9 left), we found a smooth evolution of the phase boundary Γ , while the evolution for smaller δ ≤ 0.0009 tends to be unstable. Figure 9 (right) shows a typical pattern with so called tip splitting: the evolution becomes instable and develops small fingers and side branches. For both situations we took the spatial domain to be Ω = B1 (0) the unit circle with boundary value T0 = −1, and calculated the numerical solutions using the above adaptive methods. If we choose a smaller tolerance for the adaptive method as well as finer time step sizes, the evolution is again stable, and shows qualitatively the same pattern as in the case δ = 0.006. Therefore, we believe that the observed tip splitting is due to grid effects. We now come to results obtained by simulating the anisotropic growth of a seed crystal into an undercooled melt. Starting with a circular initial interface Γ(0) := ∂B0.1 (0) should give us a feeling for the influence of the underlaying anisotropy. We calculated the evolution using the parameters δ = 0.008, A = 0.4 and k = 5. The preferred directions of this fivefold anisotropy do not coincide with the edges of the used triangulations and hence we reduce grid effects for this computation. Choosing the domain Ω = B4.0 (0) and T0 = −1 we found results as depicted in Figs.10, 11 and 12.

Fig. 11. Fivefold anisotropy: typical adaptive grid at the time t = 0.32

Using this parameters we did indeed find a transition from the initial circle to a fivefold dendrite. Figure 12 shows the graph of the temperature. Here we observe the expected behavior, nearly constant temperature inside the solid, but a fast decay in the liquid. Results of a second anisotropic calculation using the parameters δ = 0.004, A = 0.4, k = 4, T0 = −1 and the domain Ω = [−4, 4]2 are depicted in Fig. 13. According to the fourfold anisotropy the phase boundary develops four primary dendrites, we also notice very small rudiments of secondary dendrites. Choosing higher undercooling, we expect smaller temperatures inside of Ω, and such a faster growth of the solid phase. In the following experiment, we observe the evolution in case of different undercoolings. We also switched to a sixfold anisotropy, fixing their parameters to be δ = 0.01, A = 0.4 and k = 6. Figures 14 to 17 illustrate the results obtained with the boundary temperatures T0 = −1 and T0 = −2, respectively. For both examples we choose the domain Ω = B4 (0) and the initial phase boundary Γ(0) := ∂B0.1 (0).

Fig. 9. Isotropic growth for δ = 0.006 (left) and δ = 0.0009

Fig. 10. Fivefold anisotropy: Evolution of the phase boundary Γ

107

Fig. 12. Fivefold anisotropy: Graph of the temperature Θ with grid

108

M. Fried

Fig. 15. Sixfold anisotropy: Graph of the temperature with T0 = −1

Fig. 13. Anisotropic evolution of an initially circular solid phase

Figures 14 and 15 show the numerical evolution of the phase boundary Γ(t) and the graph of Θ in case of T0 = −1. As in the above examples, we observe approximately constant temperatures inside the solid and steep gradients ∇Θ in the liquid phase near the phase boundary, while the graph of Θ flattens again far away from the interface Γ(t). Reducing the boundary temperature the growth rate increases. In order to avoid an unstable evolution, we have to decrease the timestep size as well as the grid size in regions

Fig. 16. Sixfold anisotropy: Evolution of the phase boundary with T0 = −2

Fig. 17. Sixfold anisotropy: Graph of the temperature with T0 = −2

Fig. 14. Sixfold anisotropy: Evolution of the phase boundary with T0 = −1

where |∇Θ| is big. This is done by decreasing the prescribed tolerance in the above described self adaptive methods. The evolution of the solid phase obtained with the boundary temperature T0 = −2 shows some topological changes.

A level set based finite element algorithm for the simulation of dendritic growth

Fig. 18. Sixfold anisotropy: Graph of the temperature with T0 = −2, zoom to the solid phase

Fig. 19. Topological changes: Dependency on the grid size

Different fingers touch each other and merge. In between this fingers notches are developing, while the temperature gradient near the phase boundary becomes steeper. The following example illustrates the situation where two seed crystals separated by an undercooled liquid evolve and eventually merge during their growth. We found the topological change to be dependent on the the grid size of the underlying triangulation. This is demonstrated by the evolutions, which are shown in Fig. 19. The only difference between the two calculations was the prescribed tolerance for the error. While in both cases we found the evolution to be stable, prescribing a smaller tolerance leads to a situation without topological changes as we see in the right picture of Fig. 19. Both initially given crystals where of the form of (27), but the left one was rotated by 45◦ . The parameters of the anisotropy were δ = 0.006, A = 0.4 and k = 4. References 1. Bänsch, E.: Local mesh refinement in 2 and 3 dimensions, IMPACT Comput. Sci. Eng., 181–191 (1991) 2. Caginalp, G.: Length scales in phase transition models: Phase field, Cahn-Hillard and blow-up problems. In: Pattern formation: symmetry methods and applications. Chadam, J., Golubitsky, M., Langford, W., Wetton, B. (eds.), Fields Institute for Research in Mathematical Sciences, Waterloo, Canada, 1996, American Mathematical Society, Providence, Rhode Island, pp. 67–83 3. Chen, Y.-G., Giga, Y., Goto, S.: Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations. J. Differ. Geom. 33, 749–786 (1991) 4. Chopp, D.L.: Computing minimal surfaces via level set curvature flow. Journal of Computational Physics 106, 77–91 (1993)

109 5. Crandall, M.G., Lions, P.L.: Convergent difference schemes for nonlinear parabolic equations and mean curvature motion. Numer. Math. 75, 17–41 (1996) 6. Deckelnick, K.: Error analysis for a difference scheme approximating mean curvature flow, preprint. Mathematische Fakultät Freiburg, 1998 7. Deckelnick, K., Dziuk, G.: Convergence of numerical schemes for the approximation of level set solutions to the mean curvature flow. Preprint 17-00, Mathematische Fakultät Freiburg, 2000 8. Deckelnick, K., Dziuk, G.:, Numerical approximation of mean curvature flow of graphs and level sets. Preprint 04-01, Mathematische Fakultät Freiburg, 2001 9. Dörfler, W.: A convergent adaptive algorithm for Poisson’s equation. SIAM J. Numer. Anal. 33, 1106–1124 (1996) 10. Dziuk, G.: Numerical schemes for the mean curvature flow of graphs. In: Argoul, P., Fr´emond, M., Nguyen, Q.S. (eds.), IUTAM Symposium on Variations of Domains and Free-Boundary Problems in Solid Mechanics. pp. 63–70, 1999 11. Elliott, C.M.: Approximation of curvature dependent interface motion. Tech. Rep. 96/21, Centre for Mathematical Analysis and Its Applications, University of Sussex, Falmer, Brighton BN1 9QH, UK, 1996 12. Eriksson, K., Johnson, C.: Adaptive finite element methods for parabolic problems i: A linear model problem. SIAM J. Numer. Anal. 28, 43–77 (1991) 13. Evans, L.C., Spruck, J.: Motion of level sets by mean curvature I. J. Differ. Geom. 33, 635–681 (1991) 14. Evans, L.C., Spruck, J.: Motion of level sets by mean curvature II. Trans. Amer. Math. Soc. 330, 321–333 (1992) 15. Evans, L.C., Spruck, J.: Motion of level sets by mean curvature III. J. Geom. Ana. 2, 121–150 (1992) 16. Evans, L.C., Spruck, J.: Motion of level sets by mean curvature IV. J. Geom. Ana. 5, 79–116 (1995) 17. Fried, M.: Niveauflächen zur Berechnung zweidimensionaler Dendrite, PhD thesis. Institut für Angewandte Mathematik, Universität Freiburg, April 1999 18. Fried, M., Veeser, A.: Simulation and numerical analysis of dendritic growth. In: Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems. Fiedler, B. (ed.), Springer 2001 19. Gage, M.: Evolving plane curves by curvature in relative geometries. Duke Math. J. 72, 441–466 (1993) 20. Giga, Y., Goto, S., Ishii, H.: Global existence of weak solutions for interface equations coupled with diffusion equations. Siam J. Math. Anal. 23, 821–835 (1992) 21. Glicksman, M.E., Marsh, S.P.: The dendrite. In: Handbook of Crystal Growth. Hurle, D.T.J. (ed.), Vol. l, Amsterdam, London, New York, Tokyo: North Holland 1993 22. Gurtin, M.E.: Thermomechanics of evolving phase boundaries in the plane. Oxford: Clarendon Press 1993 23. Gurtin, M.E., Soner, H.M.: Some remarks on the Stefan problem with surface structure. Quarterly of Applied Mathematics 52, 291–303 (1992) 24. Huisken, G.: Non-parametric mean curvature evolution with boundary conditions. Journal of Differential Equations 77, 369-378 (1989) 25. Langer, J.S.: Instabilities and pattern formation in crystal growth. Reviews of Modern Physics 52, 1–28 (1980) 26. Lieberman, G.: The first initial-boundary value problem for quasilinear second order parabolic equations. Ann. Sci. Norm. Sup. Pisa Ser. IV., Ser. 13, 347–387 (1986) 27. Osher, S., Sethian, J.A.: Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. 79, 12–49 (1988) 28. Paolini, M.: An efficient algorithm for computing anisotropic evolution by mean curvature. In: Damlamian, A. (ed.), Proceedings of the international Conference on curvature flows and related topics held in Levico, Italy, June 27–July 2, 1994. Tokyo: Gakkotosho. GAKUTO Int. Ser., Math. Sci. Appl. 5, 1995, pp. 199–213 29. Paolini, M., Verdi, C.: Asymptotic and numerical analyses of the mean curvature flow with a space-dependent relaxation parameter. Asymptotic Anal. 5, 553–574 (1992) 30. Schmidt, A.: Die Berechnung dreidimensionaler Dendriten mit Finiten Elementen. PhD thesis, Universität Freiburg, 1993 31. Schmidt, A.: Computation of three dimensional dendrites with finite elements. Journal of Computational Physics 125, 293–312 (1996) 32. Schmidt, A., Siebert, K.G.: ALBERT – Software for scientific computations and applications. Acta Math. Univ. Comenianae 70, 105–122 (2001)

110 33. Sethian, J.A.: Level Set Methods. Cambridge University Press, 1996 34. Sethian, J.A., Strain, J.: Crystal growth and dendrite solidification. J. Comp. Phys. 98, 231–253 (1992) 35. Strain, J.: Fast tree-based redistancing for level set computations. J. Comp. Phys. 152, 648–666 (1999) 36. Sussmann, M., Smereka, P., Osher, S.J.: A level set method for computing solutions to incompressible two-phase flow. Journal of Computational Physics 114, 146–159 (1994)

M. Fried 37. Veeser, A.: Globale Existenz von Graphen unter inhomogenen Krümmungsfluß bei Dirichlet-Randbedingungen, Master’s thesis. Institut für Angewandte Mathematik, Universitat Freiburg, December 1993 38. Visintin, A.: Models of phase transitions. Vol. 28 of Progress in Nonlinear Differential Equations and Their Applications, Boston: Birkhäuser 1996 39. Walkington, N.J.: Algorithms for computing motion by mean curvature. SIAM J. Numer. Anal. 33, 2215–2238 (1996)