Adaptive finite elements with high aspect ratio for ...

Viewer
Transcript

Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576 www.elsevier.com/locate/cma

Adaptive ﬁnite elements with high aspect ratio for dendritic growth of a binary alloy including ﬂuid ﬂow induced by shrinkage Jacek Narski

*,1,

Marco Picasso

Institut d’Analyse et Calcul Scientiﬁque, Ecole Polytechnique Fe´de´rale de Lausanne, 1015 Lausanne, Switzerland Received 1 November 2005; accepted 30 October 2006 Available online 19 March 2007

Abstract An adaptive phase ﬁeld model for the solidiﬁcation of binary alloys in two space dimensions is presented. The ﬂuid ﬂow in the liquid due to diﬀerent liquid/solid densities is taken into account. The unknowns are the phase ﬁeld, the alloy concentration and the velocity/ pressure in the liquid. Continuous, piecewise linear ﬁnite elements are used for the space discretization, a semi-implicit scheme is used for the time discretization. An adaptive method allows the number of degrees of freedom to be reduced, the mesh triangles having high aspect ratio whenever needed. Numerical results are presented, showing that the ﬂuid ﬂow due to diﬀerent liquid/solid densities must be taken into account. Moreover, numerical experiments pertaining to micro-porosity and hot tearing are reported. 2007 Elsevier B.V. All rights reserved. Keywords: Mesh adaptation; Solidiﬁcation; Phase-ﬁeld; Navier-Stokes; Finite element

1. Introduction In recent years, considerable progress has been made in numerical simulation of solidiﬁcation processes at microscopic scale [1]. Although sharp interface [2,3] and levelset models [4,5] have proved to be eﬃcient, the phase ﬁeld method emerged as a method of choice in order to simulate dendritic growth in binary alloys [6–11]. In phase ﬁeld models, the location of the solid and liquid phases in the computational domain is described by introducing an order parameter, the phase ﬁeld, which varies smoothly from one in the solid to zero in the liquid through a slightly diﬀused interface. The main diﬃculty when solving numerically phase ﬁeld models is due to the very rapid change of the phase ﬁeld (and also of the concentration ﬁeld in alloys) across the diﬀuse interface, whose thickness has to be taken very small (between 1 and 10 nm) to correctly capture the *

1

Corresponding author. E-mail address: jacek.narski@epﬂ.ch (J. Narski). Supported by the Swiss National Science Foundation.

0045-7825/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.10.033

physics of the phase transformation. A high spatial resolution is therefore needed to describe the smooth transition. In order to reduce the computational time and the number of grid points adaptive isotropic ﬁnite elements [12,13] have been used. Further reduction of the number of nodes has been achieved using adaptive ﬁnite elements with high aspect ratio [14,15]. The inﬂuence of inter-dendritic liquid ﬂow has already been taken into account in dendritic simulations [16–20]. Also, the inter-dendritic liquid ﬂow induced by shrinkage – that is to say the ﬂow induced by the fact that solid and liquid densities are diﬀerent – has been considered [21–27]. The importance of extending the existing models are following. Firstly, the velocity ﬁeld may alter the solidiﬁcation process resulting in the diﬀerent shapes/sizes of the microstructures. Secondly, the pressure drop due to shrinkage, when exceeding certain critical value, can give rise to micro-porosity formation. In this paper, we address those points and study numerically the solidiﬁcation process in order to investigate the importance of the diﬀerence between solid and liquid densities. Adaptive ﬁnite elements

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

with high aspect ratio will be used as in [14,15]. Numerical results show the inﬂuence of the liquid ﬂow due to shrinkage on the formation of microstructures and demonstrate that the method is capable of predicting qualitatively micro-porosity and hot tears formation [28]. The outline of the paper is following. In the following section, we present the model. The numerical method is described in Section 3 and the numerical validation is proposed in Section 4. In Section 5, we study the inﬂuence of shrinkage on dendritic growth. While it does not alter the formation of a single solutal dendrite, it changes the interface stability of a solidiﬁcation front. The pressure drop due to shrinkage is observed during solidiﬁcation of several dendrites. The pressure drop during the hot tearing process is also investigated. 2. The model The equations governing the solidiﬁcation process are derived using a volume averaging technique in the similar way as in [16,17]. The key idea is to develop two sets of equations (for solid and liquid phases) and transform them into one set using averaging over small volume and introducing average quantities. In the following we present and discuss the averaged mass, momentum and species conservation equation for binary alloy undergoing a solid/liquid phase transition at prescribed temperature. As in [21–25], we take into account the fact that the solid and liquid densities are diﬀerent. 2.1. Mass conservation The solidiﬁcation of a binary alloy in a bounded domain X of R2 between time 0 and tend is considered. Let / : X ð0; tend Þ ! R be the phase ﬁeld describing presence of solid ð/ ¼ 1Þ or liquid ð/ ¼ 0Þ. The phase ﬁeld / varies smoothly but rapidly from zero to one in a thin region of width d, the so-called solid/liquid diﬀused interface. Let qs and ql be the constant solid and liquid densities (for most alloys ql < qs ). Then the average density q : X ð0; tend Þ ! R is deﬁned by q ¼ qs / þ ql ð1 /Þ:

ð1Þ

3563

o ðq / þ ql ð1 /ÞÞ þ divðqs /vs þ ql ð1 /Þvl Þ ¼ 0; ot s which writes, using (1) and (2): oq þ divðqvÞ ¼ 0: ot

ð3Þ

An equivalent formulation is: qs ql o/ divv ¼ þ ðv rÞ/ ot q thus the solidiﬁcation shrinkage is a source of mass in the liquid. It should be noted that, in the sharp interface limit (that is to say when the width of the solid–liquid interface d tends to zero), then the phase ﬁeld / is the characteristic function of the solid so that the density q becomes a step function and (3) has to be understood in the sense of distributions. Then, the following relation holds on the solid/ liquid interface: ½qV þ ½qv n ¼ 0; where [Æ] denotes the jump of the inside quantity across the interface, V is the normal velocity of the solid/liquid interface and the vector n denotes the normal to the interface. For instance, when solid is not moving ðvs ¼ 0Þ, this condition reduces to q ql V; ð4Þ vl n ¼ s ql see Fig. 1. At this point it should be noted that certain geometrical conﬁgurations are incompatible with the mass conservation Eq. (3). This is the case when a liquid region is surrounded by a solid region, see Fig. 2 where examples of compatible and not compatible conﬁgurations are shown. Consider the ﬁrst conﬁguration in Fig. 2. In that case the solidiﬁcation shrinkage is compensated by a liquid ﬂow through the boundary. Integrating the mass conservation Eq. (3) over the computational domain X and applying the divergence theorem gives Z Z d ðqs / þ ql ð1 /ÞÞ þ ðqs /vs þ ql ð1 /Þvl Þ n ¼ 0; dt X oX

Let vs, vl be the solid and liquid velocities, respectively. In this model, the elastic deformations of the solid are neglected. Indeed, deformation in the solid is mainly a macroscopic phenomena due to temperature variations. Therefore, the solid velocity vs is a known constant (in most cases vs ¼ 0, except for the hot tearing experiment), whereas vl is unknown. Then, the average velocity v : X ð0; tend Þ ! R2 is deﬁned by qv ¼ qs /vs þ ql ð1 /Þvl :

ð2Þ

Averaging the mass conservation equation in the solid and liquid regions yields

Fig. 1. Solidiﬁcation of a solid seed with density qs larger than the liquid density ql: as the solid/liquid interface moves with normal velocity V toward the boundary of the calculation domain X, liquid ﬂows with velocity vl toward the solid.

3564

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

Fig. 2. Examples of solidiﬁcation conﬁgurations that are compatible and not compatible with Eq. (3). The ﬁrst two conﬁgurations are compatible with (3) since some external liquid can enter the computational domain in order to feed all liquid regions. The last two conﬁgurations are not compatible with (3) since external liquid cannot feed the liquid region located at the center of the computational domain.

which becomes, since only liquid is present along the boundary oX Z d ðqs jXs ðtÞj þ ql jXl ðtÞjÞ ¼ ql vl n: ð5Þ dt oX Here above jXs ðtÞj, jXl ðtÞj are the solid and liquid area, respectively deﬁned by Z Z jXs ðtÞj ¼ /; jXl ðtÞj ¼ ð1 /Þ: X

X

Since the computational domain X is constant we obviously have jXs ðtÞj þ jXl ðtÞj ¼ X; so that (5) becomes d ql jXs ðtÞj ¼ dt qs ql

where rs and rl are the solid and liquid stress tensors, respectively. Our goal is to obtain a momentum equation for the unknown average velocity v. First, the deﬁnitions (1) and (2) together with (3) are used to eliminate the liquid velocity vl in the momentum equation: ov /qs q þ ðqv rÞv þ div q ðv vs Þ ðv vs Þ ot ð1 /Þql divð/rs þ ð1 /Þrl Þ ¼ 0: Moreover, the solid mechanical deformation is neglected, thus rs ¼ 0. Also, given a penalty parameter e 1, the penalty term 1 2 / ðv vs Þ e

Z

vl n:

oX

Thus, when solidiﬁcation occurs the ﬁrst term in the above equality is positive and, since qs > ql , then liquid must enter the computational domain to feed the liquid. Consider now the third conﬁguration in Fig. 2. Since the solid velocity vs ¼ 0, integrating the mass conservation Eq. (3) over the computational domain X now gives d ðq jXs ðtÞj þ ql jXl ðtÞjÞ ¼ 0; dt s

that is

ðqs ql Þ

d jXs ðtÞj ¼ 0: dt

Therefore, when solidiﬁcation occurs, the above equality can no longer be satisﬁed and the model is no longer valid. In practice, a third gaseous phase should appear, which leads to a micro-pore formation. Since this eﬀect is not included in the present model, the computation should be stopped when such a not compatible conﬁguration is encountered. 2.2. Momentum conservation Averaging the momentum conservation equation in the solid and liquid regions yields o ðq /vs þ ql ð1 /Þvl Þ þ divðqs /vs vs þ ql ð1 /Þvl vl Þ ot s divð/rs þ ð1 /Þrl Þ ¼ 0;

is added to the momentum equation in order to force the average velocity ﬁeld v to equal the solid velocity vs in the solid region ð/ ¼ 1Þ. This term is consistent with the method proposed in [17], where the interfacial stress term is modelled and the resulting additional term is propor2 tional to /d2 v. Finally, the liquid stress tensor ð1 /Þrl is replaced by pI þ 2ll ðvÞ where p is the average pressure, ll the liquid viscosity and ðvÞ ¼ 1=2ðrv þ rvT Þ the rate of deformation tensor of the average velocity v. Therefore, the momentum equation writes, in the whole computational domain X: ov /qs q þ ðqv rÞv þ div q ðv vs Þ ðv vs Þ ot ð1 /Þql 1 2divðll ðvÞÞ þ rp þ /2 ðv vs Þ ¼ 0: ð6Þ e It should be stressed that in the liquid far from the solid/liquid interface ð/ ¼ 0Þ, the mass and momentum Eqs. (3) and (6) reduce to the incompressible Navier–Stokes equations whereas in the solid region ð/ ¼ 1Þ, due to the penalty term, the velocity equals the solid velocity vs as e becomes small. The pressure is deﬁned in the whole computational domain but it has a physical sense only in the liquid region. In this paper, we are looking forward to predicting microporosity and hot tears formation. We now explain how these two phenomena can be related to the pressure drop along a long, thin channel, see Fig. 3. The left situation in Fig. 3 cor-

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

Cl ¼

qC kqs / þ ql ð1 /Þ

and

Cs ¼

3565

kqC : kqs / þ ql ð1 /Þ ð8Þ

Averaging the species conservation gives

Fig. 3. Pressure drop during solidiﬁcation of a liquid channel (left) and hot tearing of a liquid channel (right). In both cases, liquid can enter only from the right side of the domain.

responds to solidiﬁcation of a liquid channel, the solidiﬁcation speed of the solid–liquid interface being V. Then, from (4), we know that the velocity of the liquid at the interface is vl ¼ V ðqs ql Þ=ql , moving in the opposite direction. The right situation in Fig. 3 corresponds to a simpliﬁed experiment in which the upper and lower solid regions are (mechanically) moved with opposite directions and constant velocity vs. In both situations, we chose the boundary conditions in such a way that liquid can enter only from the right side. From mass conservation in the liquid, the liquid velocity at the right side of the domain is 2vl L=H . The pressure drop can be roughly estimated by integrating the momentum equation along the boundary of the liquid region. Assuming long channels ðL H Þ and neglecting inertial terms, the pressure drop is proportional to ll vl

L : H2

ð7Þ

Cavitation will occur whenever the pressure drop is of the order of 1 kPa. The exact value is material dependent. Therefore, if the pressure drop is suﬃcient, then microporosity should occur for the left situation of Fig. 3, whereas the right situation should predict hot tears formation, see [28,29] for the description of hot tearing. When setting ll ¼ 102 kg=ðmsÞ, H ¼ 105 m, L ¼ 103 m, then the pressure drop is 500 Pa provided vl ¼ 0:005 m=s. In the case of micro-porosity (Fig. 3, left) this would correspond to a solidiﬁcation speed of V ¼ 0:05 m=s provided the shrinkage ql =ðqs ql Þ ¼ 0:1, this being diﬃcult to observe. However, in the case of hot tearing (Fig. 3, right), this would correspond to a pulling velocity in the solid vs ’ 0:005 m=s, this being a reasonable value, see the numerical results in Section 5. 2.3. Species conservation We proceed as in [20, Eqs. (33) and (34)]. Let C be the average massic concentration of the alloy deﬁned by qC ¼ qs /C s þ ql ð1 /ÞC l ; where Cs, Cl are the solid and liquid concentrations, respectively. Introducing the constant partition coeﬃcient k k¼

Cs Cl

yields

o ðq /C s þ ql ð1 /ÞC l Þ þ divðqs /C s vs þ ql ð1 /ÞC l vl Þ ot s divðDs qs /rC s þ Dl ql ð1 /ÞrC l Þ ¼ 0; where Ds and Dl are the constant solid and liquid diﬀusion coeﬃcients. Using (8) and (2), the quantities Cs, Cl and vl can be eliminated to obtain oðqCÞ qC þ div ðqv þ ðk 1Þqs /vs Þ ot kqs / þ ql ð1 /Þ qCðql kqs Þ r/ div Dð/ÞrðqCÞ þ Dð/Þ kqs / þ ql ð1 /Þ ¼ 0; where we have set Dð/Þ ¼

kqs /Ds þ ql ð1 /ÞDl : kqs / þ ql ð1 /Þ

Introducing the volumic concentration c ¼ qC, the above equation writes oc c þ div ðqv þ ðk 1Þqs /vs Þ ot kqs / þ ql ð1 /Þ ~ /Þr/Þ divðDð/Þrc þ Dðc; ¼ 0;

ð9Þ

where we have set e /Þ ¼ Dð/Þ Dðc;

cðql kqs Þ : kqs / þ ql ð1 /Þ

2.4. Phase ﬁeld As in [17], we consider a standard phase ﬁeld equation for the solid phase moving with constant velocity vs: 1 o/ þ vs r/ lk ot /ð1 /Þð1 2/Þ ¼ C divðAðr/Þr/Þ d2 c /ð1 /Þ : ð10Þ þ T m þ ml T kqs / þ ql ð1 /Þ d Here lk denotes the kinetic mobility, C is the Gibbs–Thomson coeﬃcient. The term divAðr/Þr/ is the functional derivative (that is to say the Frechet derivative) of the surface energy Z 1 2 2 aðhðr/ðxÞÞÞ jr/ðxÞj dx; 2 X where a is the real-valued function deﬁned by

3566

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576 Table 1 Values of the physical parameters Tm

b to Fig. 4. A simple example of transformation from reference triangle K generic triangle K.

k

1000 K a

0.63 lk

0.04

0.0015 m/(K s)

Ds

Dl 10

5 10 ml

2

5 10

m =s

260 K

qs 1000 kg/m3

C 9

2

m =s ql 950 kg/m3

5 · 107 km ll 0.014 kg/(ms)

aðhÞ ¼ 1 þ a cosð4hÞ;

2.5. Summary of the model

with a the anisotropy parameter and where hðnÞ denotes the angle between a vector n 2 R2 n f0g and the ﬁrst component of the orthonormal Cartesian basis ðe1 ; e2 Þ, that is

The goal of the present model is to ﬁnd the phase ﬁeld / : X ð0; tend Þ ! R, the volumic concentration c : X ð0; tend Þ ! R, the velocity v : X ð0; tend Þ ! R2 and the pressure p : X ð0; tend Þ ! R, satisfying Eqs. (10), (9), (6) and (3). Natural boundary conditions apply on the boundary of the calculation domain X for /, c and v. Moreover, initial conditions at time t ¼ 0 must be prescribed for /, c and v. Existence and uniqueness of solutions for this model in the absence of liquid ﬂow and for suﬃciently small a (small anisotropy) are proved in [30]. A posteriori error estimation and adaptive ﬁnite elements are presented in [14]. Existence and convergence of solutions in presence of liquid ﬂow is an open problem.

cos hðnÞ ¼

n e1 : knk

Therefore the matrix AðÞ in (10) is deﬁned for n 2 R2 n f0g by aðhðnÞÞa0 ðhðnÞÞ a2 ðhðnÞÞ : AðnÞ ¼ aðhðnÞÞa0 ðhðnÞÞ a2 ðhðnÞÞ The term /ð1 /Þð1 2/Þ in (10) is the derivative of the double well which forces the phase ﬁeld to values close to zero or one. Finally, the last term in (10) is a source term accounting for the energy due to the solid–liquid phase transformation, where Tm is the melting temperature of the pure substance, ml is the slope of the liquidus line in the equilibrium phase diagram. The temperature ﬁeld T is a known linear function of space and time T ðx; tÞ ¼ T 0 þ tT_ þ Gd x;

ð11Þ

where T_ and G are given constants and d is a selected solidiﬁcation direction in space. It should be noted that the velocity ﬁeld v is not present in Eq. (10) so that, in the sharp interface limit, the classical Gibbs–Thomson relation between the interface velocity V, the curvature and the concentration c is recovered. Therefore, the eﬀect of the ﬂuid motion on Gibbs–Thomson relation is neglected. We refer to [24,25] for more general models including such eﬀects. Note that in [24] a thermal dendrite is considered and the Gibbs–Thomson law therein contains a term depending on the stress at the interface. This term however origins in the heat transfer, which is not considered in our model. Our approach is similar to [17], where the Gibbs–Thomson relation does not contain additional terms depending on pressure.

3. Numerical method Eqs. (10), (9), (6) and (3) are discretized in time using an order one semi-implicit scheme. Space discretization is based on continuous, piecewise linear ﬁnite elements on triangular adapted meshes. In order to reduce the number of degrees of freedom, the triangles may have large aspect ratio whenever needed. The reﬁnement/coarsening criterion is based on a posteriori error estimates and the adaptive algorithm has already been presented for elliptic problems [31], parabolic problems [32], Stokes problem [33], dendritic growth [14] and coalescence [15]. 3.1. The ﬁnite element method Let s be the time step, tn ¼ ns, n ¼ 0; 1; 2; . . . At each time step, given ﬁnite element approximations /n1 , cn1 , vn1 , pn1 of the phase, concentration, velocity and pressure at time tn1 , we are ﬁrst seeking for the phase /n such that

Fig. 5. Typical proﬁles of the phase ﬁeld (left), concentration ﬁeld (center) and velocity ﬁeld (right). The phase ﬁeld has values zero or one, except in the phase change region. The concentration ﬁeld changes rapidly across the phase change region, but may also vary outside. The proﬁles are taken along the diagonal of the last image in Fig. 9.

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

3567

Fig. 6. Solidiﬁcation front advancing with constant vertical velocity when d ¼ 0:25 106 m and TOL ¼ 0:15. Left: initial prescribed mesh (6332 vertexes). Middle: ﬁnal adapted mesh (5503 vertexes at the time t ¼ 1:5 s). Right: velocity ﬁeld ﬂowing toward the interface.

Table 2 Error between the sharp interface velocity and the phase ﬁeld approximation with respect to the interface width d, when TOL ¼ 0:15 d

Error

Number of vertexes

1 106 0:5 106 0:25 106

1:98 1010 1:54 1010 1:00 1010

6233 8692 10,892

Z n / /n1 w þ ðvs r/n Þw lk s X Z S 0 ð/n1 Þð/n /n1 Þ þ Sð/n1 Þ Aðr/n1 Þr/n rw þ w þC 2 d X ! Z n1 c /n1 ð1 /n1 Þ n1 T m þ ml T ð; t Þ w ¼ n1 n1 d kqs / þ ql ð1 / Þ X

ð12Þ for all test function w in the ﬁnite element space. Hereabove, Sð/Þ ¼ /ð1 /Þð1 2/Þ. The above weak formulation leads to an invertible matrix whenever 1 C 2 P 0; lk s 2d

that is s 6

2d2 ; lk C

ð13Þ

this being the classical restriction for the Allen–Cahn phase ﬁeld equations equation, see [34]. Then, we are looking for cn such that Z cn cn1 cn n n1 w þ ðk 1Þqs / vs Þ n n ðqv s X X kqs / þ ql ð1 / Þ Z Z ð14Þ rw þ D1 ð/n Þrcn rw ¼ D2 ðcn1 ; /n Þr/n rw

Z

X

X

for all test function w in the ﬁnite element space. Finally, we are seeking for the velocity vn and pressure pn such that

Fig. 8. Simpliﬁed hot tearing experiment at time t ¼ 0:01 s. Zoom of the velocity and isolines of the phase ﬁeld from 0.1 to 0.9.

Z vn vn1 þ qn ðvn1 rÞvn w qn s X Z /n qs n q ðvn1 vs Þ ðvn vs Þ : ðwÞ ð1 /n Þql X Z Z Z 1 n 2 n ð/ Þ ðv vs Þ w ¼ 0; þ ll ðvn Þ : ðwÞ pn divw þ X X X e Z X

divðqn vn Þq þ

Z

n

X

q q s

n1

X ah2 Z K qþ qn rpn rq ¼ 0 l K l tria:K

ð15Þ

ð16Þ

for all test functions w and q in the ﬁnite element space. In the above equation, we have set qn ¼ qs /n þ ql ð1 /n Þ,

Fig. 7. Simpliﬁed hot tearing experiment at time t ¼ 0:01 s. Left: adapted mesh. Middle: phase ﬁeld (blue: solid, white: liquid). Right: velocity. (For interpretation of color in this ﬁgure the reader is referred to see the web version of this article.)

3568

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

a > 0 is a dimensionless stabilisation parameter, hK denotes the size of a triangle K. When the triangle aspect ratio is moderate, then h2K can be replaced by the triangle area. When the triangle is highly stretched, hK should be the size in the minimum stretching direction, k2;K according to the notations of the next subsection, see [35] for a theoretical justiﬁcation. 3.2. Adaptive ﬁnite elements with large aspect ratio We now propose an adaptive ﬁnite element algorithm, the time step s being kept constant and such that (13)

holds. Let N be the number of time steps. The goal of the adaptive algorithm is to build successive triangulations with large aspect ratio such that the relative estimated error of the concentration c in the L2 ð0; T ; H 1 ðXÞÞ norm is close to a preset tolerance TOL. For this purpose, we introduce an error indicator which requires some further notations. This error indicator measures the error of the concentration c in the directions of maximum and minimum stretching of the triangle. The goal of the adaptive algorithm is then to equidistribute the error indicator in the directions of maximum and minimum stretching, and to align the directions of maximum and minimum stretching with the

Fig. 9. Dendritic growth of a single dendrite at times t ¼ 0:25, 0.5, 0.75 and 1 s when TOL ¼ 0:125. Left: adapted meshes. Right: concentration.

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

Fig. 10. Dendritic growth of a single dendrite. Velocity ﬁeld at time t ¼ 1 s.

directions of maximum and minimum error. We refer to [31,32,14,36,37] for theoretical justiﬁcations. b ! K be the For any triangle K of the mesh, let T K : K b aﬃne transformation which maps the reference triangle K into K. Let MK be the Jacobian of TK that is ^ þ tK : x ¼ T K ð^ xÞ ¼ M K x Since MK is invertible, it admits a singular value decomposition M K ¼ RTK KK P K , where RK and PK are orthogonal and where KK is diagonal with positive entries. In the following we set ! rT1;K 0 k1;K KK ¼ and RK ¼ ; 0 k2;K rT2;K with the choice k1;K P k2;K . A simple example of such a transformation is x1 ¼ H^x1 , x2 ¼ h^x2 , with H P h, thus

3569

Fig. 12. Dendritic growth of a single dendrite. Concentration ﬁeld at time t ¼ 1 s when TOL ¼ 0:25 (3133 vertexes), TOL ¼ 0:125 (11,726 vertexes), TOL ¼ 0:0625 (47,252 vertexes) and TOL ¼ 0:03125 (204,018 vertexes).

MK ¼

H 0 1 0 ; k1;K ¼ H ; k2;K ¼ h; r1;K ¼ ; r2;K ¼ ; 0 h 0 1

see Fig. 4. In other words r1;K and r2;K are the directions of maximum and minimum stretching, while k1;K and k2;K measure the amplitude of stretching. We now introduce cs the continuous, piecewise linear approximation in time deﬁned by cs ðx; tÞ ¼

t tn1 n tn t n1 c ðxÞ þ c ðxÞ; s s

tn1 6 t 6 tn ; x 2 X;

where cn1 , cn are computed using (14). Our simpliﬁed error indicator is then deﬁned on each time interval ½tn1 ; tn and each triangle K by

Fig. 11. Dendritic growth of a single dendrite. Zoom of the concentration, mesh and velocity at time t ¼ 1 s.

3570

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576 2

ðgn;K ðcs ÞÞ ¼

Z

tn

tn1

1 ocs ðk21;K ðrT1;K GK ðcs Þr1;K Þ 1=2 on 2 2k2;K L ðoKÞ

þ k22;K ðrT2;K GK ðcs Þr2;K ÞÞ1=2 : Here GK ðcs Þ is deﬁned by R ZZ ðg1 ðcs ÞÞ2 dx GK ðcs Þ ¼ R K gZZ ðcs ÞgZZ 2 ðcs Þdx K 1

ð17Þ

R

gZZ ðcs ÞgZZ 2 ðcs Þdx K 1 R ZZ 2 ðg2 ðcs ÞÞ dx K

!

ocs ocs where P ox and P ox are computed using the formula 1 2 1 0 P ocs 0 1 jKj ox 1 jK ocs B C tria:K P ox1 ðP Þ B P 2K C B C P 1 B C @ A¼ P B C: ocs jKj oc tria:K @ A jKj oxs P ox ðP Þ tria:K P 2K

2

n¼1

2 tria:K ðgn;K ðcs ÞÞ 2 jrcs j 0 X

RT R

6 1:25TOL:

ð18Þ

;

ox2

P 2K

PN P 0:75TOL 6

ZZ with gZZ 1 ðcs Þ, g2 ðcs Þ being the components of the so-called Zienkiewicz–Zhu error estimator 0 1 ocs ZZ ðI PÞ ox1 C g1 ðcs Þ B A; ¼ @ gZZ 2 ðcs Þ ðI PÞ ocs

2

Our adaptive algorithm aims at building triangulations with large aspect ratio such that the relative estimated error is close to a preset tolerance TOL, that is:

jK

The matrix GK ðcs Þ is an estimation of the gradient error in triangle K, therefore the term rT1;K GK ðcs Þr1;K in (17) contributes to measuring the error in the direction of the triangle’s maximum stretching and the term Z tn 1 ocs e K ðcs Þr1;K Þ1=2 k1;K ðrT1;K G 1=2 on L2 ðoKÞ tn1 2k2;K is nothing but the estimated error in the direction of the triangle’s maximum stretching.

A suﬃcient condition to satisfy (18) is to build, for each time interval ðtn1 ; tn Þ, n ¼ 1; . . . ; N , a triangulation with large aspect ratio such that Z n Z 0:752 TOL2 t 2 2 jrcs j 6 ðgn;K ðcs ÞÞ NV n tn1 X Z n Z 1:252 TOL2 t 6 jrcs j2 NV n n1 t X for all triangle K, where NVn is the number of vertexes of the mesh at time tn. We then proceed as in [31,32] to build a mesh having elements with high aspect ratio, using the BL2D mesh generator [38]. At each vertex, the estimated error in the directions of maximum and minimum stretching is equidistributed, which yields new desired values of maximum and minimum stretching. Then, the directions of maximum and minimum stretching are aligned with the directions of maximum and minimum error gradient, namely the eigenvectors of the matrix GK ðcs Þ. The reason why the mesh adaptation algorithm is based on the error estimator for the concentration ﬁeld only is the following. When computing solutal dendrites, although u, c and v vary strongly in the small region corresponding to the solid–liquid interface, the function c may also vary in other regions, whereas u does not, see Fig. 5. The adapted

Fig. 13. Dendritic growth of a single dendrite for the error tolerance TOL ¼ 0:125. Concentration ﬁeld and concentration proﬁle along the diagonal at time t ¼ 0:75 s. Left column: liquid ﬂow due to shrinkage is considered (ql =qs ¼ 0:95, 15,159 vertexes). Right column: liquid ﬂow due to shrinkage is not considered (ql =qs ¼ 1, 16,587 vertexes).

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

3571

Fig. 14. Dendritic growth of four dendrites at times t ¼ 0:025, 0.15, 0.2 and 0.225 s. Left: adapted meshes (39,442 vertexes at ﬁnal time). Middle: concentration. Right: pressure.

mesh for the concentration error estimator should also be appropriate for the phase ﬁeld and the velocity ﬁeld. Whenever needed, a simple time step adaptation algorithm can also be used. The time step s is either divided or multiplied by two if the phase ﬁeld discrepancy between time tn1 and tn exceeds imposed bounds (typically 5%).

4. Numerical validation The present model has already been validated in [14] in the case when liquid ﬂow due to shrinkage is absent, that is to say when qs ¼ ql . Our goal is therefore to validate the convection part.

3572

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

Unless speciﬁed, the values of the parameters involved in (10), (9), (6) and (3) are reported in Table 1 and correspond to an Al–Cu alloy, see [3], Table 1 column B. The temperature is given by (11) with T 0 ¼ 993:8 K, T_ ¼ 0 and G ¼ 0, the time step is s ¼ 5 104 s. In the ﬁrst test case, the mass conservation equation is validated in the sharp interface limit, more precisely, Eq. (4) is checked numerically. A horizontal solidiﬁcation front moving with constant velocity V from the bottom to the top of the computational domain is considered, see Fig. 6. The computational domain X is a square of width 0.0005 m, the front velocity is V ¼ 105 m=s. The phase ﬁeld / is imposed and given by 1 x2 Vt 0:002 1 tanh /ðx1 ; x2 ; tÞ ¼ ; 2 d the concentration, velocity and pressure being computed using the adaptive algorithm. The initial prescribed mesh, the adapted mesh and the velocity ﬁeld at ﬁnal time are shown in Fig. 6 when d ¼ 0:25 106 m and TOL ¼ 0:15. The ﬁnal time is t ¼ 1:5 s. The error (in the L2 ðXÞ norm) between the liquid velocity vl deﬁned by (4) and the corresponding phase ﬁeld approximation computed using (15) (16) is reported in Table 2. One can see that when TOL is kept ﬁxed, then the order of convergence is Oðd1=2 Þ. The goal of the second test case is to reproduce the experiment corresponding to Fig. 3, right. This test case can be seen as a simpliﬁcation of the hot tearing experiment described in [28]. The computational domain X is a rectangle with sides 0.00025 m and 0.0002 m, the channel height is H ¼ 0:00012 m, the length is L ¼ 0:00025 m. The upper

and lower solid regions are moved with opposite directions and constant velocity vs ¼ 105 m=s. The boundary conditions are such that liquid can enter only from the right side. The ﬁnal time is t ¼ 0:01 s, the time step is s ¼ 105 s, TOL ¼ 0:15. Results are reported in Figs. 7 and 8 and clearly demonstrate that hot tearing experiments can be considered with our model. The pressure drop is 0.9 Pa, which is in good agreement with the estimation (7). 5. Numerical experiments We apply this model to the growth of a single dendrite. Next we move to the case of several dendrites growing separately and observe liquid ﬂow between them and negative pressure. We compare results for diﬀerent numbers of dendrites. Then we investigate the inﬂuence of induced liquid ﬂow on interface stability. Finally we model the hot tearing experiments where two dendrites are torn apart by an external force. We realize this numerically by applying a constant velocity in one of the dendritic arms causing it to move away from the others. This artiﬁcial movement causes a signiﬁcant pressure drop between the dendritic arms which could predict the formation of hot tears. Unless explicitly stated, all the simulations are performed with the interface width d ¼ 0:5 106 m. 5.1. Single dendrite We apply our model to a single dendrite growth and compare the results when liquid ﬂow due to shrinkage is neglected ðqs ¼ ql Þ. The initial solid grain is a circle of

Fig. 15. Dendritic growth of four dendrites at time t ¼ 0:225 s. Velocity ﬁeld.

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

diameter 5 106 m of initial concentration the equilibrium liquid concentration of the phase diagram kðT T m Þ=ml . The initial concentration in the liquid phase equals 0.02 is between the liquidus concentration at equilibrium ðT T m Þ=ml ¼ 0:0238 and the solidus concentration at equilibrium kðT T m Þ=ml ¼ 0:015. The domain is a square of side 0.0002 m. The results obtained for TOL ¼ 0:125 are presented on Figs. 9–11, the maximum observed velocity being 1:2 105 m=s.

3573

In order to study the convergence of our numerical model, both d and h (or alternatively TOL) should go to zero. In practice, the interface width d is kept small and TOL varies in a certain range. Convergence with TOL is reported in Fig. 12 with d kept constant. In order to be able to perform numerical studies of the convergence of our algorithm we have multiplied the interface width d by a factor two. The obtained series of dendritic shapes seem to converge. It is impossible to decrease the error tolerance any further – the mesh with large aspect ratio for TOL ¼ 0:03125 at the ﬁnal stage has 204,018 nodes, a two times smaller tolerance would yield a mesh with approximately one million points! In Fig. 13, the shape of the dendrites are compared when liquid ﬂow due to shrinkage is neglected ðqs ¼ ql Þ or not. When liquid ﬂow is considered, the velocity ﬁeld shifts the concentration ﬁeld towards the solid region. The obtained dendrite is slightly larger when qs ¼ ql (no convection). The eﬀect of the solidiﬁcation shrinkage is small in this case. The model however will allow to investigate the pressure drop when several dendrites are involved. 5.2. Several dendrites

Fig. 16. Dendritic growth of four dendrites at time t ¼ 0:225 s. Top: zoom of the velocity. Bottom: zoom of the concentration and adapted mesh. Same scale as in Figs. 14 and 15.

The liquid ﬂow due to shrinkage is now computed around 4 or 16 dendrites. We place 2 2 ¼ 4 or 4 4 ¼ 16 dendritic seeds in the computational domain X and let the system evolve observing the liquid ﬂow and the pressure drop during the process. The distance between the seeds is always 0.0002 m so that the size of the computational domain is 0.0006 m in the 2 2 case, respectively 0.001 m in the 4 4 case. The temperature is given by (11) with T 0 ¼ 993:8 K, T_ ¼ 10 K=s and G ¼ 0. The ﬁnal time is t ¼ 0:225 s. In the mesh adaptation routine the maximum triangulation size hmax is deﬁned. The solidiﬁcation shrinkage causes the liquid to ﬂow toward the center of the square. When the dendrites are sufﬁciently big, one can observe negative pressure appearing in the almost closed central region of the computational domain X. The velocity ﬁeld is changing signiﬁcantly during the experiment, the maximum observed velocity being 0.00085 m/s. Unlike the single dendrite case, the maximum value of the velocity changes quickly. In the beginning, where distance between dendrites is big, the maximum velocity is small, namely 0.00002 m/s. However, when neighbouring dendrite tips are suﬃciently close to each other, the velocity changes rapidly in order to supply the central region of the domain with new material. When dendrites grow further, the spacing between dendrite tips of neighbouring dendrites gets smaller and that limits the amount of liquid that can pass toward the center, which causes the pressure drop. Furthermore, the concentration in the center gets much higher than in the rest of the domain. Figs. 14–16 show the evolution of the system in the case of 4 dendrites, Fig. 17 shows results for the case of 16 dendrites. The minimal values of the pressure and the maximum amplitude of the velocity are given in Table 3.

3574

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

Fig. 17. Dendritic growth of 16 dendrites at time t ¼ 0:225 s. Top left: adapted mesh (132,271 vertexes). Top right: concentration (same scale as in Fig. 14). Bottom left: velocity (same scale as in Fig. 15). Bottom right: pressure (min. 17 Pa).

Table 3 Minimum pressure and maximum velocity for dendritic growth of 4 or 16 dendrites at time t ¼ 0:2 s Dendrites

Minimum pressure (Pa)

Maximum velocity (m/s)

4 16

4.5 17

0.00015 0.00085

Table 4 Values of the physical parameters used to study interface stability Tm

k

1466 K a

0.9 lk

0.04

0.00018 m/(K s)

Ds

Dl 13

1 · 10 ml

2

m /s

210 K

C 9

5 10 qs 1000 kg/m3

2

m =s ql 920 kg/m3

2.5 · 107 km ll 0.014 kg/(ms)

5.3. Interface stability In this section we study the inﬂuence of the liquid ﬂow due to shrinkage on the stability of a planar interface. The computational domain X is a rectangle with sides 0.001 m and 0.0004 m, with the vertical solidiﬁcation front initially located at y ¼ 0:00005 m. The temperature is given by (11) with T 0 ¼ 1445 K, T_ ¼ 4 K=s and G ¼ 10000 m=s. The physical parameters involved in this experiment diﬀer from the values used in previous simulations — they correspond to an Ni–Al alloy and are given in Table 4. Initially the concentration in the solid region equals m k T I T with T I ¼ T 0 þ 0:00005 G being the initial temperaml ture at the interface. The concentration in liquid equals T I T m . The initial concentration in liquid is uniform – this ml does not present a steady state solution, where concentration in liquid decreases exponentially in the direction normal to the interface. We observe during experiments that in this particular setting the solidiﬁcation shrinkage helps

to stabilise the interface. For temperature gradient below critical value G0 sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ mT_ DC 0 G0 ¼ D the interface is not stable [39]. DC0 denotes concentration diﬀerence between liquidus and solidus at the interface. Our simulation show that this quantity depends on the ql =qs ratio. Solidiﬁcation shrinkage causes liquid ﬂow that transport material toward the interface. Hence, the solidiﬁcation starts faster when convection is present. The front advances more rapidly in the beginning phase until the steady state is reached. Therefore for a given time, the temperature at the interface is higher than in the equal density case. This results in smaller DC0 and also the value of G0. The comparison of the interfaces is shown in Fig. 18.

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

3575

Fig. 18. Interface stability. Phase ﬁeld at time t ¼ 4 s with (left) and without (right) liquid ﬂow due to shrinkage. Fig. 21. Hot tearing. Pressure ﬁeld at t ¼ 5:515 s (0.015 s after the hot tearing has started). Negative values reaching 500 Pa is observed when the upper arm is moved at velocity vs ¼ 0:005 m=s. The liquid velocity at the entrance of the channel reaches 0.04 m/s.

Fig. 19. Hot tearing. Two dendrites are torn apart by external forces (for instance due to temperature diﬀerences at macroscopic scale) causing a signiﬁcant liquid ﬂow and a pressure drop near the base of a dendrite that may lead to creation of a pore (black region).

5.4. Hot tearing Hot tearing or hot cracking is a problem commonly encountered during the casting of large freezing range alloys. The liquid ﬂow and pressure drop during the hot tearing experiment reported in [28] is now computed. Two dendrites are torn apart by an external force (see Fig. 19 for reference). We realize this numerically by applying a constant velocity in one of the dendritic arms causing it to move away from the others. This artiﬁcial movement causes a signiﬁcant pressure drop between the dendritic arms. If tearing is fast enough, the resulting pressure drop may be suﬃcient to explain hot tears formation between the pulled dendrites. In order to reproduce this process, we start the simulation as in the previous experiment. The computational domain is X ¼ 0:0012 m 0:0004 m, with the vertical solidiﬁcation front initially at y ¼ 0:00005 m. The physical parameters are as in Table 1. The planar front destabilises and dendritic growth occurs, see Fig. 20. At time t ¼ 5:5 s, when dendritic arms are long enough, we impose a vertical tearing velocity vs ¼ 0:005 m=s to the upper arm only. This

results in a maximum velocity in the liquid of 0.04 m/s near the dendritic tip. The pressure drop near the base of the dendrite reaches 500 Pa, see Fig. 21, which indicates the possibility of hot tears creation. The value obtained from the rough estimates using the Eq. (7) in Section 2.2 is 535 Pa. Thus, our model is capable of indicating the regions where micro-porosity formation might occur. The time scale of the dendritic growth and the hot tearing phenomenas are substantially diﬀerent, the latter being much shorter. Hence we use the time step adaptation routine. The adapted time step is eight times smaller than the initial one (s ¼ 5 104 s). 6. Conclusions A phase ﬁeld model for solidiﬁcation of binary alloys is presented. The liquid ﬂow due to shrinkage (diﬀerent solid and liquid densities) is considered. The model consists in solving conservation equations for the mass, momentum and concentration. The interface evolution is governed by a standard phase ﬁeld equation. Numerical simulations are performed with ﬁnite element meshes having high aspect ratio adapted in accordance with an a posteriori error estimator for the concentration equation. Numerical experiments show that the liquid ﬂow due to shrinkage has a non negligible inﬂuence. Moreover, numerical experiments are performed in order to predict micro-porosity and hot tears formation. Acknowledgements The authors wish to thank Michel Rappaz and Alain Jacot from Computational Materials Laboratory, EPFL for

Fig. 20. Hot tearing. The phase ﬁeld at times t ¼ 1 s, 2.5 s, 4 s and 5.5 s. The mesh has 22,737 vertexes at time t ¼ 5:5 s.

3576

J. Narski, M. Picasso / Comput. Methods Appl. Mech. Engrg. 196 (2007) 3562–3576

many helpful discussion concerning the physics of solidiﬁcation. References [1] W.J. Boettinger, S.R. Coriell, A.L. Greer, A. Karma, W. Kurz, M. Rappaz, R. Trivedi, Solidiﬁcation microstructures: recent developments, future directions, Acta Mater. 48 (1) (2000) 43–70. [2] D. Juric, G. Tryggvason, A front-tracking method for dendritic solidiﬁcation, J. Comput. Phys. 123 (1) (1996) 127–148. [3] A. Jacot, M. Rappaz, A pseudo-front tracking technique for the modelling of solidiﬁcation microstructures in multi-component alloys, Acta Mater. 50 (8) (2002) 1909–1926. [4] M. Fried, A level set based ﬁnite element algorithm for the simulation of dendritic growth, Comput. Vis. Sci. 7 (2) (2004) 97–110. [5] F. Gibou, R. Fedkiw, R. Caﬂisch, S. Osher, A level set approach for the numerical simulation of dendritic growth, J. Sci. Comput. 19 (1–3) (2003) 183–199. [6] J.B. Collins, H. Levine, Diﬀuse interface model of diﬀusion-limited crystal growth, Phys. Rev. B 31 (9) (1985) 6119–6122. [7] G. Caginalp, W. Xie, Phase-ﬁeld and sharp-interface alloy models, Phys. Rev. E 48 (3) (1993) 1897–1909. [8] R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Physica D 63 (1993) 410–423. [9] A. Karma, Phase-ﬁeld model of eutectic growth, Phys. Rev. E 49 (1994) 2245–2250. [10] J.A. Warren, W.J. Boettinger, Prediction of dendritic growth and microsegregation patterns in a binary alloy using the phase-ﬁeld model, Acta Metall. Mater. 43 (2) (1995) 689–703. [11] L.Q. Chen, Phase-ﬁeld models for microstructure evolutions, Annual Rev. Mater. Res. 32 (2002) 163–194. [12] A. Schmidt, Computation of three dimensional dendrites with ﬁnite elements, J. Comput. Phys. 125 (2) (1996) 293–312. [13] N. Provatas, N. Goldenfeld, J. Dantzig, Adaptive mesh reﬁnement computation of solidiﬁcation microstructures using dynamic data structures, J. Comput. Phys. 148 (1) (1999) 265–290. [14] E. Burman, M. Picasso, Anisotropic, adaptive ﬁnite elements for the computation of a solutal dendrite, Interfaces Free Bound 5 (2) (2003) 103–127. [15] E. Burman, A. Jacot, M. Picasso, Adaptive ﬁnite elements with high aspect ratio for the computation of coalescence using a phase-ﬁeld model, J. Comput. Phys. 195 (1) (2004) 153–174. [16] J. Ni, C. Beckermann, A volume-averaged 2-phase model for transport phenomena during solidiﬁcation, Metall. Trans. B 22 (3) (1991) 349–361. [17] C. Beckermann, H.J. Diepers, I. Steinbach, A. Karma, X. Tong, Modeling melt convection in phase-ﬁeld simulations of solidiﬁcation, J. Comput. Phys. 154 (2) (1999) 468–496. [18] B. Nestler, A.A. Wheeler, L. Ratke, C. Sto¨cker, Phase-ﬁeld model for solidiﬁcation of a monotectic alloy with convection, Physica D 141 (1–2) (2000) 133–154. [19] J.H. Jeong, N. Goldenfeld, J. Dantzig, Phase ﬁeld model for threedimensional dendritic growth with ﬂuid ﬂow, Phys. Rev. E 64 (4) (2001) 041602.

[20] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Phase-ﬁeld simulation of solidiﬁcation, Annu. Rev. Mater. Res. 32 (2002) 163– 194. [21] D.M. Anderson, G.B. McFadden, A.A. Wheeler, A phase-ﬁeld model of solidiﬁcation with convection, Physica D 135 (1–2) (2000) 175–194. [22] D.M. Anderson, G.B. McFadden, A.A. Wheeler, A phase-ﬁeld model with convection: sharp-interface asymptotics, Physica D 151 (2–4) (2001) 305–331. [23] J.C. Heinrich, D.R. Poirier, Convection modeling in directional solidiﬁcation, C.R. Mecanique 332 (2004) 429–445. [24] M. Griebel, W. Merz, T. Neunhoeﬀer, Mathematical modeling and numerical simulation of freezing processes of a supercooled melt under consideration of density changes, Comp. Vis. Sci. 1 (4) (1999) 201–219. [25] Y. Sun, C. Beckermann, Diﬀuse interface modeling of two-phase ﬂows based on averaging: mass and momentum equations, Physica D 198 (3–4) (2004) 281–308. [26] M. Conti, Density change eﬀects on crystal growth from the melt, Phys. Rev. E 64 (5) (2001) 051601. [27] M. Conti, Advection ﬂow eﬀects in the growth of a free dendrite, Phys. Rev. E (Statist. Nonlinear, Soft Matter Phys.) 69 (2) (2004) 022601. [28] I. Farup, J.M. Drezet, M. Rappaz, In situ observation of hat-tearing formation in succinonitrile–acetome, Acta Mater. 49 (2) (2001) 1261– 1269. [29] M. Rappaz, J.M. Drezet, M. Gremaud, A new hot-tearing criterion, Metall. Mater. Trans. A 30 (2) (1999) 449–455. [30] E. Burman, J. Rappaz, Existence of solutions to an anisotropic phaseﬁeld model, Math. Methods Appl. Sci. 26 (13) (2003) 1137–1160. [31] M. Picasso, Numerical study of the eﬀectivity index for an anisotropic error indicator based on Zienkiewicz–Zhu error estimator, Comm. Numer. Methods Engrg. 19 (1) (2003) 13–23. [32] M. Picasso, An anisotropic error indicator based on Zienkiewicz–Zhu error estimator: application to elliptic and parabolic problems, SIAM J. Sci. Comput. 24 (4) (2003) 1328–1355, electronic. [33] M. Picasso, An adaptive algorithm for the stokes problem using continuous, piecewise linear stabilized ﬁnite elements and meshes with high aspect ratio, Appl. Numer. Math. 54 (3–4) (2005) 470–490. [34] C.M. Elliott, Approximation of curvature dependent interface motion, in: The State of the Art in Numerical Analysis (York, 1996), Oxford University Press, New York, 1997, pp. 407–440. [35] S. Micheletti, S. Perotto, M. Picasso, Stabilized ﬁnite elements on anisotropic meshes: a priori error estimates for the advection– diﬀusion and the Stokes problems, SIAM J. Numer. Anal. 41 (3) (2003) 1131–1162, electronic. [36] L. Formaggia, S. Perotto, New anisotropic a priori error estimates, Numer. Math. 89 (4) (2001) 641–667. [37] L. Formaggia, S. Perotto, Anisotropic error estimates for elliptic problems, Numer. Math. 94 (2003) 67–92. [38] H. Borouchaki, P. Laug, The bl2d mesh generator: beginner’s guide, user’s and programmer’s manua l. Technical Report RT-0194, Institut National de Recherche en Informatique et Automatique (INRIA), Rocquencourt, 78153 Le Chesnay, France, 1996. [39] W. Kurtz, D.J. Fisher, Fundamentals of Solidiﬁcation, Trans Tech Publications, Uetikon-Zuerich, 1998.