- Email: [email protected]
Phase field simulations of stress distributions in solidification structures
ARTICLE IN PRESS
Journal of Crystal Growth 310 (2008) 1331–1336 www.elsevier.com/locate/jcrysgro
Phase field simulations of stress distributions in solidification structures Takuya Ueharaa,, Motoshi Fukuib,1, Nobutada Ohnoc a
Department of Mechanical Systems Engineering, Yamagata University, 4-3-16, Jonan, Yonezawa, Yamagata 992-8510, Japan b Nagoya University, Japan c Department of Computational Science and Engineering, Nagoya University, Furo-Cho, Chikusa-ku, Nagoya 464-8603, Japan Available online 23 December 2007
Abstract Phase field simulations of stress evolution during solidification processes are carried out, resulting in the generation of complicated stress distributions in solidification microstructures. The governing equations, which include the coupling relations among phase transformation, temperature and stress/strain, are numerically solved using the two-dimensional finite element method. To avoid complexity, mechanical properties of the liquid are applied under a simplified assumption that sufficiently small elastic constants can approximate the liquid behavior. A dendritic microstructure is considered in this study as a typical solidification microstructure. As a result, a complicated stress distribution in a dendrite is obtained. Especially, strong stresses are revealed to distribute at the bottom of sidebranches of a dendrite. Subsequently, a microstructure constructed with two dendrites is simulated, and it reveals that high stresses are generated in the regions where the liquid remains till the very last stage of the solidification. Finally, two types of dendritic patterns are investigated, and the stresses are revealed to strongly depend on the morphology of the microstructures. r 2007 Elsevier B.V. All rights reserved. PACS: 81.10.Jt; 62.20.Fe; 81.40.Lm; 46.15.x Keywords: A1. Stresses; A1. Computer simulation; A1. Phase field model; A1. Solidification; A1. Dendrites
1. Introduction The microstructure of a material affects the macroscopic mechanical properties of the material. Especially in highly developed devices, such as micromechanical and electronic devices, a fine feature has a relatively large effect. For precise and reliable designing, therefore, a proper estimation of the effects of microstructures on mechanical behavior is required. Microstructures are usually formed during solidification processes and several defects, which are presumed to be caused by the stress distribution, are generated. The measurement of the stresses in microstructure is, however, difficult to be achieved, and the precise stress distributions have not been clarified so far. Computer simulation is, therefore, the only way for predicting them. The phase field model [1] is one of the most Corresponding author. Tel./fax: +81 238 26 3285.
E-mail addresses: [email protected] (T. Uehara), [email protected] (N. Ohno). 1 Present affiliation: Denso Corp. 0022-0248/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2007.12.035
promising tools for such purposes. This model has been successfully applied for demonstrating the complicated morphology of solidified materials such as dendrites in a pure material [2], dendrites in an alloy system [3], directional solidification patterns [4], facets [5], and polycrystals [6–8]. Since the fundamental equations are derived in a thermodynamically consistent way [9], including the anisotropic detail [10], this model is recognized as a physically proper model. For application to the mechanical behavior, Uehara et al. calculated elastic stresses of solid during solidification process by a one-dimensional model [11]. In order to estimate more precise stresses, they formulated the coupling effects among phase, temperature and stress/strain including plastic stress [12,13], and executed numerical simulations on the stress evolution during an in-solid phase transformation [13,14]. However, complicated solidification structures such as dendrites were not considered there. In this paper, therefore, the previous model is extended to the solidification problem, and the stress distributions in dendritic microstructures are investigated.
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1332
Mechanical equilibrium:
2. Fundamental equations In the phase field model, the present state in phase at a position is represented by a scalar phase-field variable f, and the spatial distribution and temporal evolution are expressed by a differential equation. In a two-phase problem, for example, f is defined as 0 in liquid and 1 in solid, and smoothly changes in a very thin interfacial region. In this paper, the fundamental equations derived by Uehara et al. [13] are applied. These equations are derived based on thermodynamic relations, as originally formulated for the coupling of phase and thermal fields by Wang et al. [1,9]. Mechanical coupling is included by considering elastic strain energy and dissipation due to plastic deformation. The equations used in this paper are the phase field equation, heat conduction equation, and stress equilibrium equation with elasto-plastic constitutive equations. In all relations, specifically in anisotropy and mechanical constitutive relations, two-dimensional formulas are presented. Note that the anisotropy for crystal growth is taken into account to realize the dendritic patterns, while mechanically isotropic properties are assumed. Phase field equation: _ ¼ ar ðRðyÞrfÞ þ gfð1 fÞðf 1 þ MÞ, mf 2 where RðyÞ ¼
r2 ðyÞ
rðyÞr0 ðyÞ
rðyÞr0 ðyÞ
r2 ðyÞ
(1)
! ,
(2)
rðyÞ ¼ 1 þ dr cosð4ðy y0 ÞÞ,
(3)
M ¼ bfð1 fÞðLð1 T=T m Þ þ f ðsij ÞÞ.
(4)
Here, m, a, g and b are phase-field parameters depending on the interfacial thickness, interfacial energy, etc., y is an angle expressing the normal direction of the interface, y0 is the reference direction, dr is the anisotropy strength, L is the latent heat, and T m is the temperature of equilibrium transformation. Originally, the anisotropy is introduced in the left-hand side of Eq. (1) as that for the kinetic coefficient [10], but m is assumed to be constant in this paper. A function f ðsij Þ expresses the stress dependency of a phase transformation, which is neglected in this paper for simplicity. Heat conduction equation: _ Tas_ ij dij þ sij e_ p . rcT_ ¼ kr2 T þ Lp0 ðfÞf ij
(5)
Here, r, c and k are the mass density, specific heat and heat conductivity, respectively. The function p0 ðfÞ is a differential of pðfÞ with respect to f, where pðfÞ is a function expressing the properties in the interfacial region. In this paper, pðfÞ ¼ f3 ð10 15f þ 6f2 Þ
(6)
is used. Another formula for this function is discussed in Ref. [9].
qsij ¼ 0. qxj
(7)
Here, the stress sij is assumed to be represented by a conventional elasto-plastic continuum mechanical framework. The total strain eij is expressed as eij ¼ eeij þ epij , eeij
(8) epij
where and are the elastic and plastic strains, respectively. The properties of a liquid are assumed to be approximated by taking the elastic modulus to be sufficiently small. The elastic strain is then expressed as n¯ 1 þ n¯ _ ij . _ ij þ bp0 ðfÞfd (9) s_ ij s_ kk dij þ a¯ Td ¯ E¯ E ¯ n¯ and a¯ are Young’s modulus, Poisson’s ratio and Here, E, the thermal expansion coefficient, respectively, of the model, which are given as e_eij ¼
1 1 1 ¼ ð1 pðfÞÞ þ pðfÞ , El Es E¯ ð1 pðfÞÞnl =E l þ pðfÞns =E s n¯ ¼ , ð1 pðfÞÞ=E l þ pðfÞ=E s a¯ ¼ ð1 pðfÞÞal þ pðfÞas ,
ð10Þ ð11Þ ð12Þ
where the subscripts s and l refer to values for solid and liquid, respectively. Another material parameter b in Eq. (9) is the dilatation transformation coefficient. The description of plastic strain in a solid has been discussed from a conventional continuum mechanical standpoint, and many constitutive equations and yielding conditions have been suggested [15]. Since a detailed investigation on the choice of the models is not of great importance here, an often-used model based on the normal flow law, Prager’s kinematic hardening rule, and Mises’ yielding function, presented below, is applied in this paper: e_pij ¼
9 ðsmn amn Þs_ mn ðsij aij Þ. H0 4s¯ 20
(13)
Here, sij , aij , s¯ 0 and H 0 are the deviatoric stress, deviatoric back stress, initial yield stress, and hardening coefficient, respectively, and the equivalent stress, seq ¼ ð3sij sij =2Þ1=2
(14)
is used for quantitative stress evaluation. 3. Model and conditions The governing equations described in the previous section are solved numerically using the finite element method. A square region in two dimensions is divided into 4-node isoparametric elements. The total numbers of nodes and elements are 6400 and 6561, respectively, and the mechanical boundary conditions are illustrated in Fig. 1. In the phase field model, remeshing techniques [16] are often used to improve the computational efficiency, but the
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
initially divided mesh is used throughout the simulation to avoid complication in the coding. A solid nucleus is set in a supercooled liquid as the initial state. In order to take the initial volumetric change due to the nucleation into account, the phase transformation from liquid to solid is employed using the initial 10 time steps, without solving the phase field equation. For the thermal equation, an adiabatic condition is imposed during the crystal growth process. Since the temperature in the model increases due to the latent heat generation during solidification, the model is cooled down after proper time steps with heat convection from the boundary so that solidification is complete in the whole model. For the phase field, the Neumann condition is applied: the gradient of f is set to be zero at the boundaries.
1333
The parameters applied in this study is summarized in Table 1, where T 0 is the initial temperature, h and T w are the heat convection coefficient and environmental temperature for the duration that the convectional boundary is applied, Dx is the length of an element, and Dt is the time increment. The negative value of the transformation coefficient b represents the contraction due to solidification. The relationship between the parameters and the experimental properties have been investigated previously [17–19]. This set of parameters, however, relates to a type of virtual material, and a specific material is not considered. 4. Results and discussion 4.1. Dendrite formation
y x Fig. 1. Simulation region for the FEM analysis.
First, a single dendrite formation process is simulated. A nucleus is set on the bottom-left corner, and the preferable direction of the crystal growth, which corresponds to y0 in Eq. (3), is set to 0 (and 90 due to four-fold symmetry). Note that the boundaries along the x and y axes correspond to line-symmetric boundaries under this condition. Figs. 2(a)–(c) represent the phase field, temperature, and equivalent stress distribution as time progresses. The solid area starts growing from the nucleus along the boundary. The shape of the solid region is circular at the very early stage, before the time step t ¼ 100. The latent heat
Table 1 Applied parameters for simulations m, J s=m3 r, kg=m3 T m, K E s , GPa a, 1/K Dx, m
6:93 102 1:00 103 1728 200 2:5 106 1:0 108
a, J/m c, J/(kg K) T 0, K ns b Dt, s
3:02 108 5:42 103 1511h, W=ðm2 KÞ 0.3 s¯ 0 , MPa 1:0 1011
g, J=m3 k, W/(m K) 1:0 104 E l , Pa 250 dr
3:27 108 8:50 101 T w, K 1.0 H 0 , MPa 0.015
b, m3 =J L, J=m3 1511 nl 250
Fig. 2. Variation of (a) phase field, (b) temperature and (c) equivalent stress during dendritic growth.
9:18 108 2:35 109 0.5
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1334
generation due to the phase transformation increases the temperature at the interface, and an instability is induced there. A dendritic pattern is then generated as shown in Figs. 2(a)(ii)–(vi). Thermal dilatation due to latent heat and contraction due to solidification cause the complicated stress distribution. Especially at the interfacial region, a severe condition is imposed by their almost simultaneous occurrence, and the stress value exceeds the yielding stress. When side branches are generated, high stresses are observed at the bottom of the branches, as shown in Figs. 2(c)(iv)–(vi). In these parts, solidification is completed earlier than in the adjacent area, and the stressed state has been in a plastic state. When solidification occurs in the vicinity, contraction occurs while the thermal gradient becomes more severe. Consequently, high stresses are induced in the region. 4.2. Residual stress generation When a dendrite is formed at t ¼ 1500, the temperature rises as shown in Fig. 2(vi). In order to urge solidification in the whole model, heat is taken from the boundary. A convectional heat condition is imposed on the two edges of x ¼ 0 and y ¼ 0 (left and bottom edges). The phase field, temperature, and equivalent stress are shown in Figs. 3(a)–(c), respectively. As shown in figures (a) and (b), the temperature decreases from the cooling edges, from where solidification is driven. As the phase transformation progresses, the primary branches thicken, and finally the transformation is complete in the whole region. At the beginning of the cooling, two primary branches along the x and y boundaries reach the right and top walls, respectively. This results in high stresses as shown in Fig. 3(c)(i). As the primary branches become thicker, the stress distribution varies. High stresses are concentrated near the surface at the bottom of the side branches. These regions are major solidification fronts in
the cooling process, and they are where large contractions due to the phase transformation occur, generating high stress. When the solidification is complete in the whole region, a residual stress distribution remains as shown in Fig. 3(c)(vi). High stresses are observed in the middle of the model, corresponding to the area where the solidification occurred belatedly. The top and left boundaries can be regarded as the collision of the solid regions, due to the comparability of the Neumann condition and line symmetry, which are discussed in the next section. 4.3. Two-dendrite models In the previous section, a single dendrite was considered. In a realistic situation, a microstructure is composed of multiple dendrites. The influence of the collision of two dendrites on the stress distribution is discussed in this section. Two nuclei are set at the bottom left and top right corners. The preferable growth direction is set as y0 ¼ 0 for both dendrites, while another combination is applied later. The simulation region is set larger than before, with 10 000 nodes and 10 201 elements. The other conditions are the same as in the previous section, except that all boundaries are set as cooling edges in the cooling stage. The stress evolution in the growing process is almost identical to that in Fig. 2, because the collision of the two dendrites does not occur. The stress in each dendrite does not affect the other since they are separated by the liquid. The collision occurs during the cooling process, hence only the results after t ¼ 1500 are presented. Figs. 4(a) and (b) represent the results of the phase field and equivalent stress, respectively. As time progresses, the remaining liquid between the two dendrites gradually solidifies. The liquid region shrinks inward, as shown in Figs. 4(ii)–(v), and the transformation is finally complete at the center. The residual stress distribution is exhibited in Fig. 4(b)(vi).
Fig. 3. Variation of (a) phase field, (b) temperature and (c) equivalent stress during cooling process.
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1335
Fig. 4. Variation of (a) phase field and (b) equivalent stress during cooling process for mono-orientated dendrites.
Fig. 5. Variation of (a) phase field and (b) equivalent stress during cooling process for offset orientated dendrites.
The high stress is distributed along the diagonal lines. One of the diagonals, from the bottom left to the top right, corresponds to the 45 zone caused by the four-fold symmetry of each dendrite. The other diagonal, from the top left to the bottom right, is along the liquid region remaining between the two dendrites as in Fig. 4(i). Both these regions are where the liquid remained in the later stage; it is concluded that high stresses are generated where the solidification occurs belatedly compared to the adjacent regions. How the microstructural pattern affects the stress distribution is now investigated. The reference orientation y0 is set as 0 for one nucleus and 45 for a second nucleus. All other conditions are identical to those in the former simulation. Fig. 5 represents the result of the phase field and equivalent stresses. Similar to Fig. 4, a complicated stress distribution in the solidified area is observed, and large residual stresses are generated in the model. The residual stresses in Figs. 4(vi) and 5(vi) are drawn in Fig. 6 with a magnified range so that the distribution can be seen more clearly. While the phase is uniform for both cases, the residual stress distribution is obviously different. Also the maximum values of the stress for (a) and (b) are different in that they are 260 and 295 MPa, respectively. Therefore, it is concluded that the residual stresses are strongly affected by the morphology of the microstructure.
Fig. 6. Comparison of the residual stress distributions for two types of microstructures. (a) Same orientation; (b) offset orientation.
5. Concluding remarks The stress evolution and residual stress distribution in a solidified material are simulated using the phase field model, and the following results were obtained:
(1) Complicated stresses are generated in dendritic microstructures. (2) Residual stresses are retained in solidified material. (3) Stress distribution is strongly dependent on the morphology of the microstructure. Since these stresses are difficult to be measured experimentally, the model and procedure described in this paper will be effective for estimating the stresses in the solidified materials.
ARTICLE IN PRESS 1336
T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
There are still many challenges for improving the present model. For instance, the flow and viscosity of liquid is not taken into account. Mechanical anisotropy needs to be considered with respect to the anisotropy in crystal growth from a crystallographic viewpoint. A conventional plastic constitutive relation is assumed in this paper. Essentially, the source of the plasticity, such as dynamics of dislocations and sliding behavior, should be inclusively introduced. Three-dimensional modeling is also necessary for quantitative evaluation, as well as the use of realistic parameters. As future works, a physical understanding of the mechanical behavior both in elasticity and plasticity especially at the interfaces, and its formulation in the phase field framework are required. The authors have started an investigation based on molecular dynamics simulations [20], which is believed to be an appropriate approach to achieve a comprehensive multi-scale modeling of the plastic behavior of materials. References [1] R.F. Sekerka, Fundamentals of phase field theory, in: Advances in Crystal Growth Research, Elsevier, Amsterdam, 2001, p. 21. [2] R. Kobayashi, Physica D 63 (1993) 410.
[3] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys. Rev. E 47 (1993) 1893. [4] Z. Bi, R.F. Sekerka, Physica A 261 (1998) 95. [5] T. Uehara, R.F. Sekerka, J. Crystal Growth 254 (2003) 251. [6] D. Fan, L.-Q. Chen, Acta Mater. 45 (1997) 611. [7] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415. [8] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Acta Mater. 51 (2003) 6035. [9] S.-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, G.B. McFadden, Physica D 69 (1993) 189. [10] G.B. McFadden, A.A. Wheeler, R.J. Braun, S.R. Coriell, J. Crystal Growth 48 (1993) 2016. [11] T. Uehara, T. Tsujino, J. Crystal Growth 275 (2005) e219. [12] T. Uehara, T. Tsujino, Trans. JSME A 72 (2006) 438. [13] T. Uehara, T. Tsujino, N. Ohno, J. Crystal Growth 300 (2007) 530. [14] T. Uehara, T. Tsujino, Trans. JSME A 72 (2006) 849. [15] For example: R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, London, 1951. [16] N. Provatas, N. Goldenfeld, J. Dantzig, J. Comp. Phys. 148 (1999) 265. [17] G. Caginalp, P. Fife, Phys. Rev. B 33 (1986) 7792. [18] B.T. Murray, A.A. Wheeler, M.E. Glicksman, J. Crystal Growth 154 (1995) 386. [19] A. Karma, W.J. Rappel, Phys. Rev. E 57 (1998) 4323. [20] T. Uehara, N. Wakabayashi, N. Ohno, Key Eng. Mater. 340/341 (2007) 1003.
Journal of Crystal Growth 310 (2008) 1331–1336 www.elsevier.com/locate/jcrysgro
Phase field simulations of stress distributions in solidification structures Takuya Ueharaa,, Motoshi Fukuib,1, Nobutada Ohnoc a
Department of Mechanical Systems Engineering, Yamagata University, 4-3-16, Jonan, Yonezawa, Yamagata 992-8510, Japan b Nagoya University, Japan c Department of Computational Science and Engineering, Nagoya University, Furo-Cho, Chikusa-ku, Nagoya 464-8603, Japan Available online 23 December 2007
Abstract Phase field simulations of stress evolution during solidification processes are carried out, resulting in the generation of complicated stress distributions in solidification microstructures. The governing equations, which include the coupling relations among phase transformation, temperature and stress/strain, are numerically solved using the two-dimensional finite element method. To avoid complexity, mechanical properties of the liquid are applied under a simplified assumption that sufficiently small elastic constants can approximate the liquid behavior. A dendritic microstructure is considered in this study as a typical solidification microstructure. As a result, a complicated stress distribution in a dendrite is obtained. Especially, strong stresses are revealed to distribute at the bottom of sidebranches of a dendrite. Subsequently, a microstructure constructed with two dendrites is simulated, and it reveals that high stresses are generated in the regions where the liquid remains till the very last stage of the solidification. Finally, two types of dendritic patterns are investigated, and the stresses are revealed to strongly depend on the morphology of the microstructures. r 2007 Elsevier B.V. All rights reserved. PACS: 81.10.Jt; 62.20.Fe; 81.40.Lm; 46.15.x Keywords: A1. Stresses; A1. Computer simulation; A1. Phase field model; A1. Solidification; A1. Dendrites
1. Introduction The microstructure of a material affects the macroscopic mechanical properties of the material. Especially in highly developed devices, such as micromechanical and electronic devices, a fine feature has a relatively large effect. For precise and reliable designing, therefore, a proper estimation of the effects of microstructures on mechanical behavior is required. Microstructures are usually formed during solidification processes and several defects, which are presumed to be caused by the stress distribution, are generated. The measurement of the stresses in microstructure is, however, difficult to be achieved, and the precise stress distributions have not been clarified so far. Computer simulation is, therefore, the only way for predicting them. The phase field model [1] is one of the most Corresponding author. Tel./fax: +81 238 26 3285.
E-mail addresses: [email protected] (T. Uehara), [email protected] (N. Ohno). 1 Present affiliation: Denso Corp. 0022-0248/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2007.12.035
promising tools for such purposes. This model has been successfully applied for demonstrating the complicated morphology of solidified materials such as dendrites in a pure material [2], dendrites in an alloy system [3], directional solidification patterns [4], facets [5], and polycrystals [6–8]. Since the fundamental equations are derived in a thermodynamically consistent way [9], including the anisotropic detail [10], this model is recognized as a physically proper model. For application to the mechanical behavior, Uehara et al. calculated elastic stresses of solid during solidification process by a one-dimensional model [11]. In order to estimate more precise stresses, they formulated the coupling effects among phase, temperature and stress/strain including plastic stress [12,13], and executed numerical simulations on the stress evolution during an in-solid phase transformation [13,14]. However, complicated solidification structures such as dendrites were not considered there. In this paper, therefore, the previous model is extended to the solidification problem, and the stress distributions in dendritic microstructures are investigated.
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1332
Mechanical equilibrium:
2. Fundamental equations In the phase field model, the present state in phase at a position is represented by a scalar phase-field variable f, and the spatial distribution and temporal evolution are expressed by a differential equation. In a two-phase problem, for example, f is defined as 0 in liquid and 1 in solid, and smoothly changes in a very thin interfacial region. In this paper, the fundamental equations derived by Uehara et al. [13] are applied. These equations are derived based on thermodynamic relations, as originally formulated for the coupling of phase and thermal fields by Wang et al. [1,9]. Mechanical coupling is included by considering elastic strain energy and dissipation due to plastic deformation. The equations used in this paper are the phase field equation, heat conduction equation, and stress equilibrium equation with elasto-plastic constitutive equations. In all relations, specifically in anisotropy and mechanical constitutive relations, two-dimensional formulas are presented. Note that the anisotropy for crystal growth is taken into account to realize the dendritic patterns, while mechanically isotropic properties are assumed. Phase field equation: _ ¼ ar ðRðyÞrfÞ þ gfð1 fÞðf 1 þ MÞ, mf 2 where RðyÞ ¼
r2 ðyÞ
rðyÞr0 ðyÞ
rðyÞr0 ðyÞ
r2 ðyÞ
(1)
! ,
(2)
rðyÞ ¼ 1 þ dr cosð4ðy y0 ÞÞ,
(3)
M ¼ bfð1 fÞðLð1 T=T m Þ þ f ðsij ÞÞ.
(4)
Here, m, a, g and b are phase-field parameters depending on the interfacial thickness, interfacial energy, etc., y is an angle expressing the normal direction of the interface, y0 is the reference direction, dr is the anisotropy strength, L is the latent heat, and T m is the temperature of equilibrium transformation. Originally, the anisotropy is introduced in the left-hand side of Eq. (1) as that for the kinetic coefficient [10], but m is assumed to be constant in this paper. A function f ðsij Þ expresses the stress dependency of a phase transformation, which is neglected in this paper for simplicity. Heat conduction equation: _ Tas_ ij dij þ sij e_ p . rcT_ ¼ kr2 T þ Lp0 ðfÞf ij
(5)
Here, r, c and k are the mass density, specific heat and heat conductivity, respectively. The function p0 ðfÞ is a differential of pðfÞ with respect to f, where pðfÞ is a function expressing the properties in the interfacial region. In this paper, pðfÞ ¼ f3 ð10 15f þ 6f2 Þ
(6)
is used. Another formula for this function is discussed in Ref. [9].
qsij ¼ 0. qxj
(7)
Here, the stress sij is assumed to be represented by a conventional elasto-plastic continuum mechanical framework. The total strain eij is expressed as eij ¼ eeij þ epij , eeij
(8) epij
where and are the elastic and plastic strains, respectively. The properties of a liquid are assumed to be approximated by taking the elastic modulus to be sufficiently small. The elastic strain is then expressed as n¯ 1 þ n¯ _ ij . _ ij þ bp0 ðfÞfd (9) s_ ij s_ kk dij þ a¯ Td ¯ E¯ E ¯ n¯ and a¯ are Young’s modulus, Poisson’s ratio and Here, E, the thermal expansion coefficient, respectively, of the model, which are given as e_eij ¼
1 1 1 ¼ ð1 pðfÞÞ þ pðfÞ , El Es E¯ ð1 pðfÞÞnl =E l þ pðfÞns =E s n¯ ¼ , ð1 pðfÞÞ=E l þ pðfÞ=E s a¯ ¼ ð1 pðfÞÞal þ pðfÞas ,
ð10Þ ð11Þ ð12Þ
where the subscripts s and l refer to values for solid and liquid, respectively. Another material parameter b in Eq. (9) is the dilatation transformation coefficient. The description of plastic strain in a solid has been discussed from a conventional continuum mechanical standpoint, and many constitutive equations and yielding conditions have been suggested [15]. Since a detailed investigation on the choice of the models is not of great importance here, an often-used model based on the normal flow law, Prager’s kinematic hardening rule, and Mises’ yielding function, presented below, is applied in this paper: e_pij ¼
9 ðsmn amn Þs_ mn ðsij aij Þ. H0 4s¯ 20
(13)
Here, sij , aij , s¯ 0 and H 0 are the deviatoric stress, deviatoric back stress, initial yield stress, and hardening coefficient, respectively, and the equivalent stress, seq ¼ ð3sij sij =2Þ1=2
(14)
is used for quantitative stress evaluation. 3. Model and conditions The governing equations described in the previous section are solved numerically using the finite element method. A square region in two dimensions is divided into 4-node isoparametric elements. The total numbers of nodes and elements are 6400 and 6561, respectively, and the mechanical boundary conditions are illustrated in Fig. 1. In the phase field model, remeshing techniques [16] are often used to improve the computational efficiency, but the
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
initially divided mesh is used throughout the simulation to avoid complication in the coding. A solid nucleus is set in a supercooled liquid as the initial state. In order to take the initial volumetric change due to the nucleation into account, the phase transformation from liquid to solid is employed using the initial 10 time steps, without solving the phase field equation. For the thermal equation, an adiabatic condition is imposed during the crystal growth process. Since the temperature in the model increases due to the latent heat generation during solidification, the model is cooled down after proper time steps with heat convection from the boundary so that solidification is complete in the whole model. For the phase field, the Neumann condition is applied: the gradient of f is set to be zero at the boundaries.
1333
The parameters applied in this study is summarized in Table 1, where T 0 is the initial temperature, h and T w are the heat convection coefficient and environmental temperature for the duration that the convectional boundary is applied, Dx is the length of an element, and Dt is the time increment. The negative value of the transformation coefficient b represents the contraction due to solidification. The relationship between the parameters and the experimental properties have been investigated previously [17–19]. This set of parameters, however, relates to a type of virtual material, and a specific material is not considered. 4. Results and discussion 4.1. Dendrite formation
y x Fig. 1. Simulation region for the FEM analysis.
First, a single dendrite formation process is simulated. A nucleus is set on the bottom-left corner, and the preferable direction of the crystal growth, which corresponds to y0 in Eq. (3), is set to 0 (and 90 due to four-fold symmetry). Note that the boundaries along the x and y axes correspond to line-symmetric boundaries under this condition. Figs. 2(a)–(c) represent the phase field, temperature, and equivalent stress distribution as time progresses. The solid area starts growing from the nucleus along the boundary. The shape of the solid region is circular at the very early stage, before the time step t ¼ 100. The latent heat
Table 1 Applied parameters for simulations m, J s=m3 r, kg=m3 T m, K E s , GPa a, 1/K Dx, m
6:93 102 1:00 103 1728 200 2:5 106 1:0 108
a, J/m c, J/(kg K) T 0, K ns b Dt, s
3:02 108 5:42 103 1511h, W=ðm2 KÞ 0.3 s¯ 0 , MPa 1:0 1011
g, J=m3 k, W/(m K) 1:0 104 E l , Pa 250 dr
3:27 108 8:50 101 T w, K 1.0 H 0 , MPa 0.015
b, m3 =J L, J=m3 1511 nl 250
Fig. 2. Variation of (a) phase field, (b) temperature and (c) equivalent stress during dendritic growth.
9:18 108 2:35 109 0.5
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1334
generation due to the phase transformation increases the temperature at the interface, and an instability is induced there. A dendritic pattern is then generated as shown in Figs. 2(a)(ii)–(vi). Thermal dilatation due to latent heat and contraction due to solidification cause the complicated stress distribution. Especially at the interfacial region, a severe condition is imposed by their almost simultaneous occurrence, and the stress value exceeds the yielding stress. When side branches are generated, high stresses are observed at the bottom of the branches, as shown in Figs. 2(c)(iv)–(vi). In these parts, solidification is completed earlier than in the adjacent area, and the stressed state has been in a plastic state. When solidification occurs in the vicinity, contraction occurs while the thermal gradient becomes more severe. Consequently, high stresses are induced in the region. 4.2. Residual stress generation When a dendrite is formed at t ¼ 1500, the temperature rises as shown in Fig. 2(vi). In order to urge solidification in the whole model, heat is taken from the boundary. A convectional heat condition is imposed on the two edges of x ¼ 0 and y ¼ 0 (left and bottom edges). The phase field, temperature, and equivalent stress are shown in Figs. 3(a)–(c), respectively. As shown in figures (a) and (b), the temperature decreases from the cooling edges, from where solidification is driven. As the phase transformation progresses, the primary branches thicken, and finally the transformation is complete in the whole region. At the beginning of the cooling, two primary branches along the x and y boundaries reach the right and top walls, respectively. This results in high stresses as shown in Fig. 3(c)(i). As the primary branches become thicker, the stress distribution varies. High stresses are concentrated near the surface at the bottom of the side branches. These regions are major solidification fronts in
the cooling process, and they are where large contractions due to the phase transformation occur, generating high stress. When the solidification is complete in the whole region, a residual stress distribution remains as shown in Fig. 3(c)(vi). High stresses are observed in the middle of the model, corresponding to the area where the solidification occurred belatedly. The top and left boundaries can be regarded as the collision of the solid regions, due to the comparability of the Neumann condition and line symmetry, which are discussed in the next section. 4.3. Two-dendrite models In the previous section, a single dendrite was considered. In a realistic situation, a microstructure is composed of multiple dendrites. The influence of the collision of two dendrites on the stress distribution is discussed in this section. Two nuclei are set at the bottom left and top right corners. The preferable growth direction is set as y0 ¼ 0 for both dendrites, while another combination is applied later. The simulation region is set larger than before, with 10 000 nodes and 10 201 elements. The other conditions are the same as in the previous section, except that all boundaries are set as cooling edges in the cooling stage. The stress evolution in the growing process is almost identical to that in Fig. 2, because the collision of the two dendrites does not occur. The stress in each dendrite does not affect the other since they are separated by the liquid. The collision occurs during the cooling process, hence only the results after t ¼ 1500 are presented. Figs. 4(a) and (b) represent the results of the phase field and equivalent stress, respectively. As time progresses, the remaining liquid between the two dendrites gradually solidifies. The liquid region shrinks inward, as shown in Figs. 4(ii)–(v), and the transformation is finally complete at the center. The residual stress distribution is exhibited in Fig. 4(b)(vi).
Fig. 3. Variation of (a) phase field, (b) temperature and (c) equivalent stress during cooling process.
ARTICLE IN PRESS T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
1335
Fig. 4. Variation of (a) phase field and (b) equivalent stress during cooling process for mono-orientated dendrites.
Fig. 5. Variation of (a) phase field and (b) equivalent stress during cooling process for offset orientated dendrites.
The high stress is distributed along the diagonal lines. One of the diagonals, from the bottom left to the top right, corresponds to the 45 zone caused by the four-fold symmetry of each dendrite. The other diagonal, from the top left to the bottom right, is along the liquid region remaining between the two dendrites as in Fig. 4(i). Both these regions are where the liquid remained in the later stage; it is concluded that high stresses are generated where the solidification occurs belatedly compared to the adjacent regions. How the microstructural pattern affects the stress distribution is now investigated. The reference orientation y0 is set as 0 for one nucleus and 45 for a second nucleus. All other conditions are identical to those in the former simulation. Fig. 5 represents the result of the phase field and equivalent stresses. Similar to Fig. 4, a complicated stress distribution in the solidified area is observed, and large residual stresses are generated in the model. The residual stresses in Figs. 4(vi) and 5(vi) are drawn in Fig. 6 with a magnified range so that the distribution can be seen more clearly. While the phase is uniform for both cases, the residual stress distribution is obviously different. Also the maximum values of the stress for (a) and (b) are different in that they are 260 and 295 MPa, respectively. Therefore, it is concluded that the residual stresses are strongly affected by the morphology of the microstructure.
Fig. 6. Comparison of the residual stress distributions for two types of microstructures. (a) Same orientation; (b) offset orientation.
5. Concluding remarks The stress evolution and residual stress distribution in a solidified material are simulated using the phase field model, and the following results were obtained:
(1) Complicated stresses are generated in dendritic microstructures. (2) Residual stresses are retained in solidified material. (3) Stress distribution is strongly dependent on the morphology of the microstructure. Since these stresses are difficult to be measured experimentally, the model and procedure described in this paper will be effective for estimating the stresses in the solidified materials.
ARTICLE IN PRESS 1336
T. Uehara et al. / Journal of Crystal Growth 310 (2008) 1331–1336
There are still many challenges for improving the present model. For instance, the flow and viscosity of liquid is not taken into account. Mechanical anisotropy needs to be considered with respect to the anisotropy in crystal growth from a crystallographic viewpoint. A conventional plastic constitutive relation is assumed in this paper. Essentially, the source of the plasticity, such as dynamics of dislocations and sliding behavior, should be inclusively introduced. Three-dimensional modeling is also necessary for quantitative evaluation, as well as the use of realistic parameters. As future works, a physical understanding of the mechanical behavior both in elasticity and plasticity especially at the interfaces, and its formulation in the phase field framework are required. The authors have started an investigation based on molecular dynamics simulations [20], which is believed to be an appropriate approach to achieve a comprehensive multi-scale modeling of the plastic behavior of materials. References [1] R.F. Sekerka, Fundamentals of phase field theory, in: Advances in Crystal Growth Research, Elsevier, Amsterdam, 2001, p. 21. [2] R. Kobayashi, Physica D 63 (1993) 410.
[3] A.A. Wheeler, W.J. Boettinger, G.B. McFadden, Phys. Rev. E 47 (1993) 1893. [4] Z. Bi, R.F. Sekerka, Physica A 261 (1998) 95. [5] T. Uehara, R.F. Sekerka, J. Crystal Growth 254 (2003) 251. [6] D. Fan, L.-Q. Chen, Acta Mater. 45 (1997) 611. [7] R. Kobayashi, J.A. Warren, W.C. Carter, Physica D 119 (1998) 415. [8] J.A. Warren, R. Kobayashi, A.E. Lobkovsky, W.C. Carter, Acta Mater. 51 (2003) 6035. [9] S.-L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, G.B. McFadden, Physica D 69 (1993) 189. [10] G.B. McFadden, A.A. Wheeler, R.J. Braun, S.R. Coriell, J. Crystal Growth 48 (1993) 2016. [11] T. Uehara, T. Tsujino, J. Crystal Growth 275 (2005) e219. [12] T. Uehara, T. Tsujino, Trans. JSME A 72 (2006) 438. [13] T. Uehara, T. Tsujino, N. Ohno, J. Crystal Growth 300 (2007) 530. [14] T. Uehara, T. Tsujino, Trans. JSME A 72 (2006) 849. [15] For example: R. Hill, The Mathematical Theory of Plasticity, Oxford University Press, London, 1951. [16] N. Provatas, N. Goldenfeld, J. Dantzig, J. Comp. Phys. 148 (1999) 265. [17] G. Caginalp, P. Fife, Phys. Rev. B 33 (1986) 7792. [18] B.T. Murray, A.A. Wheeler, M.E. Glicksman, J. Crystal Growth 154 (1995) 386. [19] A. Karma, W.J. Rappel, Phys. Rev. E 57 (1998) 4323. [20] T. Uehara, N. Wakabayashi, N. Ohno, Key Eng. Mater. 340/341 (2007) 1003.