Adaptive finite elements with high aspect ratio for dendritic growth of a binary alloy including fluid flow induced by shrinkage

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

Adaptive finite elements with high aspect ratio for dendritic growth of a binary alloy including fluid flow induced by shrinkage Jacek Narski

*,1,

Marco Picasso

Institut d’Analyse et Calcul Scientifique, 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 field model for the solidification of binary alloys in two space dimensions is presented. The fluid flow in the liquid due to different liquid/solid densities is taken into account. The unknowns are the phase field, the alloy concentration and the velocity/ pressure in the liquid. Continuous, piecewise linear finite 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 fluid flow due to different 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; Solidification; Phase-field; Navier-Stokes; Finite element

1. Introduction In recent years, considerable progress has been made in numerical simulation of solidification processes at microscopic scale [1]. Although sharp interface [2,3] and levelset models [4,5] have proved to be efficient, the phase field method emerged as a method of choice in order to simulate dendritic growth in binary alloys [6–11]. In phase field models, the location of the solid and liquid phases in the computational domain is described by introducing an order parameter, the phase field, which varies smoothly from one in the solid to zero in the liquid through a slightly diffused interface. The main difficulty when solving numerically phase field models is due to the very rapid change of the phase field (and also of the concentration field in alloys) across the diffuse 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@epfl.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 finite elements [12,13] have been used. Further reduction of the number of nodes has been achieved using adaptive finite elements with high aspect ratio [14,15]. The influence of inter-dendritic liquid flow has already been taken into account in dendritic simulations [16–20]. Also, the inter-dendritic liquid flow induced by shrinkage – that is to say the flow induced by the fact that solid and liquid densities are different – has been considered [21–27]. The importance of extending the existing models are following. Firstly, the velocity field may alter the solidification process resulting in the different 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 solidification process in order to investigate the importance of the difference between solid and liquid densities. Adaptive finite 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 influence of the liquid flow 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 influence of shrinkage on dendritic growth. While it does not alter the formation of a single solutal dendrite, it changes the interface stability of a solidification front. The pressure drop due to shrinkage is observed during solidification of several dendrites. The pressure drop during the hot tearing process is also investigated. 2. The model The equations governing the solidification 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 different. 2.1. Mass conservation The solidification 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 field describing presence of solid ð/ ¼ 1Þ or liquid ð/ ¼ 0Þ. The phase field / varies smoothly but rapidly from zero to one in a thin region of width d, the so-called solid/liquid diffused 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 defined 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 solidification 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 field / 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 configurations 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 configurations are shown. Consider the first configuration in Fig. 2. In that case the solidification shrinkage is compensated by a liquid flow 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 defined by qv ¼ qs /vs þ ql ð1  /Þvl :

ð2Þ

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

Fig. 1. Solidification 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 flows with velocity vl toward the solid.

3564

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

Fig. 2. Examples of solidification configurations that are compatible and not compatible with Eq. (3). The first two configurations are compatible with (3) since some external liquid can enter the computational domain in order to feed all liquid regions. The last two configurations 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 defined 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 definitions (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 solidification occurs the first 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 configuration 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 solidification occurs, the above equality can no longer be satisfied and the model is no longer valid. In practice, a third gaseous phase should appear, which leads to a micro-pore formation. Since this effect is not included in the present model, the computation should be stopped when such a not compatible configuration 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 field 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 defined 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 solidification 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 solidification of a liquid channel, the solidification 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 simplified 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 sufficient, 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 solidification speed of V ¼ 0:05 m=s provided the shrinkage ql =ðqs  ql Þ ¼ 0:1, this being difficult 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 defined by qC ¼ qs /C s þ ql ð1  /ÞC l ; where Cs, Cl are the solid and liquid concentrations, respectively. Introducing the constant partition coefficient 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 diffusion coefficients. 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 field As in [17], we consider a standard phase field 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 coefficient. 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 defined 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 first component of the orthonormal Cartesian basis ðe1 ; e2 Þ, that is

The goal of the present model is to find the phase field / : 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 flow and for sufficiently small  a (small anisotropy) are proved in [30]. A posteriori error estimation and adaptive finite elements are presented in [14]. Existence and convergence of solutions in presence of liquid flow is an open problem.

cos hðnÞ ¼

n  e1 : knk

Therefore the matrix AðÞ in (10) is defined 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 field 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 field 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 solidification direction in space. It should be noted that the velocity field 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 effect of the fluid motion on Gibbs–Thomson relation is neglected. We refer to [24,25] for more general models including such effects. 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 finite 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 refinement/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 finite element method Let s be the time step, tn ¼ ns, n ¼ 0; 1; 2; . . . At each time step, given finite element approximations /n1 , cn1 , vn1 , pn1 of the phase, concentration, velocity and pressure at time tn1 , we are first seeking for the phase /n such that

Fig. 5. Typical profiles of the phase field (left), concentration field (center) and velocity field (right). The phase field has values zero or one, except in the phase change region. The concentration field changes rapidly across the phase change region, but may also vary outside. The profiles 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. Solidification front advancing with constant vertical velocity when d ¼ 0:25  106 m and TOL ¼ 0:15. Left: initial prescribed mesh (6332 vertexes). Middle: final adapted mesh (5503 vertexes at the time t ¼ 1:5 s). Right: velocity field flowing toward the interface.

Table 2 Error between the sharp interface velocity and the phase field 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 finite 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 field 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 finite element space. Finally, we are seeking for the velocity vn and pressure pn such that

Fig. 8. Simplified hot tearing experiment at time t ¼ 0:01 s. Zoom of the velocity and isolines of the phase field 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 finite element space. In the above equation, we have set qn ¼ qs /n þ ql ð1  /n Þ,

Fig. 7. Simplified hot tearing experiment at time t ¼ 0:01 s. Left: adapted mesh. Middle: phase field (blue: solid, white: liquid). Right: velocity. (For interpretation of color in this figure 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 justification. 3.2. Adaptive finite elements with large aspect ratio We now propose an adaptive finite 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 field at time t ¼ 1 s.

directions of maximum and minimum error. We refer to [31,32,14,36,37] for theoretical justifications. b ! K be the For any triangle K of the mesh, let T K : K b affine 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 field 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 defined 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 simplified error indicator is then defined 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 defined 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 sufficient 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 field 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 field and concentration profile along the diagonal at time t ¼ 0:75 s. Left column: liquid flow due to shrinkage is considered (ql =qs ¼ 0:95, 15,159 vertexes). Right column: liquid flow 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 final time). Middle: concentration. Right: pressure.

mesh for the concentration error estimator should also be appropriate for the phase field and the velocity field. 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 field 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 flow 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 specified, 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 first test case, the mass conservation equation is validated in the sharp interface limit, more precisely, Eq. (4) is checked numerically. A horizontal solidification 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 field / 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 field at final time are shown in Fig. 6 when d ¼ 0:25 106 m and TOL ¼ 0:15. The final time is t ¼ 1:5 s. The error (in the L2 ðXÞ norm) between the liquid velocity vl defined by (4) and the corresponding phase field approximation computed using (15) (16) is reported in Table 2. One can see that when TOL is kept fixed, 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 simplification 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 final 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 flow between them and negative pressure. We compare results for different numbers of dendrites. Then we investigate the influence of induced liquid flow 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 artificial movement causes a significant 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 flow 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 field.

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 final 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 flow due to shrinkage is neglected ðqs ¼ ql Þ or not. When liquid flow is considered, the velocity field shifts the concentration field towards the solid region. The obtained dendrite is slightly larger when qs ¼ ql (no convection). The effect of the solidification 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 flow 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 flow 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 final time is t ¼ 0:225 s. In the mesh adaptation routine the maximum triangulation size hmax is defined. The solidification shrinkage causes the liquid to flow toward the center of the square. When the dendrites are sufficiently big, one can observe negative pressure appearing in the almost closed central region of the computational domain X. The velocity field is changing significantly 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 sufficiently 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 influence of the liquid flow 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 solidification 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 differ 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 solidification shrinkage helps

to stabilise the interface. For temperature gradient below critical value G0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mT_ DC 0 G0 ¼ D the interface is not stable [39]. DC0 denotes concentration difference between liquidus and solidus at the interface. Our simulation show that this quantity depends on the ql =qs ratio. Solidification shrinkage causes liquid flow that transport material toward the interface. Hence, the solidification 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 field at time t ¼ 4 s with (left) and without (right) liquid flow due to shrinkage. Fig. 21. Hot tearing. Pressure field 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 differences at macroscopic scale) causing a significant liquid flow 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 flow 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 artificial movement causes a significant pressure drop between the dendritic arms. If tearing is fast enough, the resulting pressure drop may be sufficient 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 solidification 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 different, 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 field model for solidification of binary alloys is presented. The liquid flow due to shrinkage (different 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 field equation. Numerical simulations are performed with finite 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 flow due to shrinkage has a non negligible influence. 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 field 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 solidification. References [1] W.J. Boettinger, S.R. Coriell, A.L. Greer, A. Karma, W. Kurz, M. Rappaz, R. Trivedi, Solidification microstructures: recent developments, future directions, Acta Mater. 48 (1) (2000) 43–70. [2] D. Juric, G. Tryggvason, A front-tracking method for dendritic solidification, J. Comput. Phys. 123 (1) (1996) 127–148. [3] A. Jacot, M. Rappaz, A pseudo-front tracking technique for the modelling of solidification microstructures in multi-component alloys, Acta Mater. 50 (8) (2002) 1909–1926. [4] M. Fried, A level set based finite element algorithm for the simulation of dendritic growth, Comput. Vis. Sci. 7 (2) (2004) 97–110. [5] F. Gibou, R. Fedkiw, R. Caflisch, 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, Diffuse interface model of diffusion-limited crystal growth, Phys. Rev. B 31 (9) (1985) 6119–6122. [7] G. Caginalp, W. Xie, Phase-field 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-field 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-field model, Acta Metall. Mater. 43 (2) (1995) 689–703. [11] L.Q. Chen, Phase-field models for microstructure evolutions, Annual Rev. Mater. Res. 32 (2002) 163–194. [12] A. Schmidt, Computation of three dimensional dendrites with finite elements, J. Comput. Phys. 125 (2) (1996) 293–312. [13] N. Provatas, N. Goldenfeld, J. Dantzig, Adaptive mesh refinement computation of solidification microstructures using dynamic data structures, J. Comput. Phys. 148 (1) (1999) 265–290. [14] E. Burman, M. Picasso, Anisotropic, adaptive finite elements for the computation of a solutal dendrite, Interfaces Free Bound 5 (2) (2003) 103–127. [15] E. Burman, A. Jacot, M. Picasso, Adaptive finite elements with high aspect ratio for the computation of coalescence using a phase-field model, J. Comput. Phys. 195 (1) (2004) 153–174. [16] J. Ni, C. Beckermann, A volume-averaged 2-phase model for transport phenomena during solidification, 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-field simulations of solidification, J. Comput. Phys. 154 (2) (1999) 468–496. [18] B. Nestler, A.A. Wheeler, L. Ratke, C. Sto¨cker, Phase-field model for solidification of a monotectic alloy with convection, Physica D 141 (1–2) (2000) 133–154. [19] J.H. Jeong, N. Goldenfeld, J. Dantzig, Phase field model for threedimensional dendritic growth with fluid flow, Phys. Rev. E 64 (4) (2001) 041602.

[20] W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Phase-field simulation of solidification, Annu. Rev. Mater. Res. 32 (2002) 163– 194. [21] D.M. Anderson, G.B. McFadden, A.A. Wheeler, A phase-field model of solidification with convection, Physica D 135 (1–2) (2000) 175–194. [22] D.M. Anderson, G.B. McFadden, A.A. Wheeler, A phase-field 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 solidification, C.R. Mecanique 332 (2004) 429–445. [24] M. Griebel, W. Merz, T. Neunhoeffer, 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, Diffuse interface modeling of two-phase flows based on averaging: mass and momentum equations, Physica D 198 (3–4) (2004) 281–308. [26] M. Conti, Density change effects on crystal growth from the melt, Phys. Rev. E 64 (5) (2001) 051601. [27] M. Conti, Advection flow effects 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 phasefield model, Math. Methods Appl. Sci. 26 (13) (2003) 1137–1160. [31] M. Picasso, Numerical study of the effectivity 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 finite 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 finite elements on anisotropic meshes: a priori error estimates for the advection– diffusion 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 Solidification, Trans Tech Publications, Uetikon-Zuerich, 1998.