1 Introduction

Dendritic crystal structures which form during solidification of metallic alloys have their enormous influence on the mechanical properties of solid alloys, thus, the dendritic growth problem has been a topic of long-term interest within the academia and metal industry. An exact solution for the non-dimensional temperature distribution ahead of an isolated dendritic crystal, having the form of an isothermal and semi-infinite paraboloid of revolution with a fixed radius of tip curvature, ρ, and growing at a constant velocity, V, was suggested by Papapetrou in 1935 and first presented by Ivantsov [1] in 1947.

Later Horvay and Cahn [2] provided a rigorous mathematical solution of the steady-state diffusion process around parabolic interfaces including the 3D paraboloid of revolution, the elliptical paraboloid, and the 2D parabolic plate. The ability to correctly predict ρ is a problem of fundamental importance to the theory of dendritic growth with an additional dendritic tip selection constant (stability parameter) σ* needed to determine the operating conditions (the combination of radius ρ and growth rate V) at the dendrite tip. Values of the σ* have been calculated based on two main dendritic growth theories: marginal stability [3] and microscopic solvability [4]. The marginal stability estimates the tip selection parameter as constant for all materials under all conditions to be 1/4π2 ≈ 0.025 which is very close to the numerical value of σ* = 0.025 evaluated by Langer and Müller-Krumbhaar [3] for the symmetric problem in 3D, but approximately twice smaller than a value calculated using 3D phase-field simulations by Oguchi and Suzuki [5] for Al-4.5wt%Cu as a one-sided problem (the solute diffusion in the solid phase is negligible). Based on the marginal stability theory the Kurz-Giovanola-Trivedi (KGT) constrained (columnar) dendritic growth model was developed for steady-state conditions with the stability parameter which is equal to 0.0253 [6]. It must be emphasized that the KGT model or similar models, like the LKT (Lipton, Kurz and Trivedi) one, are an inherent part of a range of alloy solidification models like cellular automaton (CA) [7], front tracking method (FTM) [8] or one-domain multiphase models based on volume averaging [9]. Gäumann et al. [10] have numerically modified Hunt’s columnar to equiaxed transition (CET) analytical model [11] by combining the KGT model. Rebow and Browne [12] estimated dendritic tip stability parameters σ* of two aluminium alloys, namely Al-4wt%Cu and Al-2wt%Si, based on measured values of their crystal-melt surface energy anisotropy strength ε and a simple linear scaling law of microscopic solvability theory. Subsequently they showed that the stability parameter has the significant influence on columnar dendritic growth models and consequently on a columnar to equiaxed transition with help of modified Hunt’s analytical maps and meso-scale front tracking simulations. Recently, Mullis [13] presented results from a phase field model and found that the tip stability parameter σ* is non constant, but varies as a function of tip undercooling for different alloy concentration, Lewis number, and equilibrium partition coefficient.

Usually the tip stability parameter σ*, which is experimentally determined for 3D situations, is applied to 2D models. Altundas and Caginalp [14] showed that the 2D and 3D phase field simulations of the tip velocity differ by a factor of approximately 1.9.

In the presented study, authors attempt to cure the situation by deriving 2D dendritic tip scaling parameter for aluminium-based alloy Al-4wt%Cu. The obtained parameter is then incorporated into the KGT dendritic growth model in order to compare it with the original 3D KGT counterpart and to derive two-dimensional and three-dimensional versions of the modified Hunt’s analytical model for the columnar-to-equiaxed transition (CET). The conclusions drawn from the above analysis are further confirmed through numerical calculations of the two cases of a Al-4wt%Cu metallic alloy solidification using the front tracking technique on a fixed control-volume grid. Results, including the mush-liquid front position, the calculated solutal under-cooling, the average temperature gradient at a front of the dendrite tip envelope and a new predictor of the relative tendency to form an equiaxed zone, are shown, compared and discussed.

2 Dendritic growth model and columnar to Equiaxed transition (CET) analytical map

Following the same approach developed by Rebow and Browne [12] dendritic tip stability parameters σ* of Al-4wt%Cu for 2D equal to 0.02 and 3D equal 0.058 based on a simple linear scaling law of microscopic solvability theory, were derived. These values of σ*, along with the typical 3D value of 0.0253 for the symmetric problem used in dendritic growth models were introduced into the KGT model with diffusion fields in 2D and 3D described by 2 different functions, namely 3D Ivantsov and 2D Horvay & Cahn and with the following properties used for Al-4wt%Cu: diffusivity of solute in liquid (D l) 2.4 × 10−9 m2/s; the Gibbs-Thomson coefficient (Γ) 2.36 × 10−7 mK; the partition coefficient (k) 0.17; the liquidus slope (m) -2.6 K/ wt%. The dendrite tip velocity V vs. the solutal (constitutional) tip undercooling ΔT c for Al-4wt%Cu is presented in Fig. 1 and fitted to the following relationship: V = A (ΔT c )n where A ms−1 K-n and n are fitting parameters, which values are shown in that figure. The dendrite tip velocity V for the 2D dendritic tip stability parameter is much lower than for 3D ones.

Fig. 1
figure 1

Comparison of the KGT dendritic growth model for 2D and 3D stability parameters σ* of Al-4wt%Cu alloy

Full size image

A Columnar to Equiaxed Transition (CET) analytical map can be significantly changed due to a selection of dendritic growth models as proved by Rebow and Browne [12]. The sensitivity of CET maps for 2D and 3D different stability parameters incorporated into the KGT model is compared. Based on a modification of Hunt’s model [11], a relationship between the temperature gradient, G and the dendrite tip solutal undercooling, ΔT c for fully equiaxed growth for 2D and 3D can be written as,

$$ \begin{array}{l}{G}_{2 D}={M}_{2 D}\;{N}_0^{1/2}\Delta {T}_c\kern0.5em \left(1-\frac{\Delta {T}_N^{n+1}}{\Delta {T}_c^{n+1}}\right)\\ {}{G}_{3 D}={M}_{3 D}\;{N}_0^{1/3}\Delta {T}_c\kern0.5em \left(1-\frac{\Delta {T}_N^{n+1}}{\Delta {T}_c^{n+1}}\right)\end{array} $$
(1)

where N 0 is the total number of heterogeneous nucleation sites per unit volume or area, \( {N}_{0\left(3 D\right)}\left[1/{m}^3\right]=\frac{\pi}{6}{N}_{0\left(2 D\right)}^{3/2}\left[1/{m}^2\right] \), ΔT N is nucleation undercooling, M is constant for the particular growth model:

$$ \begin{array}{l}{M}_{2 D}=\frac{1}{n+1}\sqrt{\frac{-\pi}{ \ln \left(1-\Phi \right)}}\\ {}{M}_{3 D}=\frac{1}{n+1}\sqrt[3]{\frac{-4\pi}{3 \ln \left(1-\Phi \right)}}\end{array} $$
(2)

where n is a parameter, see eq. (1) and the volume fraction of equiaxed grains ϕ = 0.49 for a 3D map [11] and ϕ = 0.544 for a 2D map.

Substituting V = A (ΔT c )n into eq. (1) along with values of the parameter n, ΔT N  = 0.75 K, N 02D = 1.24 106 m−2 and N 03D = 109 m−3, the CET map presented in Fig. 2 for Al-4wt%Cu alloy and 2D and 3D stability parameters are established.

Fig. 2
figure 2

CET map for an Al-4wt%Cu alloy

Full size image

In general, results show that for the same G value, the CET shifts to a higher V value for the 3D σ* stability parameter proving that there is a significant difference in 2D and 3D representation of the CET effect.

3 Mathematical and numerical models

The volume averaged heat balance equation, describing the transient, diffusive heat transfer in the solidifying binary alloy is discretised on the control volume mesh. The integral on the single control volume (CV) reads

$$ \underset{V}{\int}\frac{\partial \left(\rho cT\right)}{\partial t} dV=\underset{\partial V}{\int } k\frac{\partial T}{\partial {x}_j}{n}_j dS+\underset{V}{\int}\rho L\frac{\partial {g}_S}{\partial t} dV $$
(3)

where ρ, c, k and L are density, specific heat, thermal conductivity and latent heat, respectively. All material properties are assumed constant and equal in both phases. The last term of the right side of the eq. (3) is due to heat release accompanying phase change. The solute micro-diffusion model obeys the Scheil’s model, so the volumetric solid fraction g S can be expressed as a function of temperature T:

$$ {g}_S=\left\{\begin{array}{l}0\operatorname{}\operatorname{}\operatorname{} for\kern0.5em {T}_{liq}\le T\\ {}1-{\left(\frac{T_M- T}{T_M-{T}_{liq}}\right)}^{\frac{1}{k_p-1}}\operatorname{} for\kern0.5em {T}_E\le T<{T}_{liq}\\ {}1\operatorname{}\operatorname{}\operatorname{} for\kern0.5em T<{T}_E\end{array}\right. $$
(4)

where T M , T E , T liq , k p are the melting temperature of the solvent, the eutectic temperature, the liquidus temperature and the partition coefficient.

Between the liquidus and solidus isotherms, a region of coexistence of solid and liquid phases develops. Typically, two dendrite morphologies are observed, namely motionless columnar dendrites fixed to cooled walls of the mould and equiaxed grains in the central part of the domain. To distinguish these regions the virtual surface, representing an envelope of columnar dendrites tips is defined which is moving across the domain with respect to predefined columnar dendrite tips growth kinetics. Such defined interface is represented with mass-less markers connected with line segments (Fig. 3). Their motion mimics the growth of columnar dendrites tips and is dependent on the local under-cooling ∆T, which is a difference between the liquidus temperature and the interpolated temperature at the markers. The bi-linear interpolation scheme is used to determine temperatures in front markers on the basis of nodal temperatures got from adjacent cells centres (black dots in Fig. 3). Positions of the markers are known from the previous time step, they are denoted as black dots connected with black segments in Fig. 3. The new positions of markers (red dots in the Fig. 3) shifted along the normal vectors to front, are calculated with the formula

$$ {\mathbf{X}}_i^n={\mathbf{X}}_{i,0}+{\mathbf{n}}_i\cdot {w}_i\cdot \Delta t $$
(5)
Fig. 3
figure 3

A fragment of the front crossing the control volume mesh in previous (black line) and next time step (red line).

Full size image

where X i , 0 is the position of the i-th marker in the previous time step, \( {\mathbf{X}}_i^n \) is the position of the moved i-th marker in the n-th iteration, n i is the vector normal to the front determined in the i-th marker, w i dendrite growth rate calculated in the i-th marker according to the KGT model, and ∆t is the time step. Details of the procedure of the front tracking are given in [8, 15].

The originally proposed by Browne [16] predictor for equiaxed solidification is determined to investigate the impact of the dendrite tip kinetics on the tendency to formation of CET. In the author’s opinion the blockage of growth of columnar dendrites by thickening equiaxed grains growing in the under-cooled liquid could take place till the equiaxed index parameter achieves its maximum value. The definition of equiaxed index is taken from [16]

$$ EI(t)=\sum_i{U}_b(i)\Delta V, $$
(6)

where U b is the level of under-cooling,

$$ {U}_b(i)=\left\{\begin{array}{l}0\operatorname{}\operatorname{} if\kern0.5em T(i)>{T}_{liq}\\ {}{T}_{liq}- T(i)\operatorname{} if\kern0.5em T(i)\le {T}_{liq}\kern0.5em and\kern0.5em I(i)<0.5\\ {}0\operatorname{}\operatorname{} if\kern0.5em T(i)\le {T}_{liq}\kern0.5em and\kern0.5em I(i)\ge 0.5\end{array}\right., $$
(7)

related to temperature T(i) and the indicator function I(i) in the i-th control volume. The latter is equal to 1 behind the front and 0 in the rest of the domain, namely in front of the dendrite tips envelope. The immersed front tracking approach, originally introduced by Peskin [17] is utilized and the procedure for determination of the I(x) function was adopted after [17, 18]. The step variation across the interface is relaxed using the Pesking distribution functions in the vicinity of the interface. The resulting grid-gradient field G(x) which has nonzero values only close to the front, is the right hand side of the Poisson equation

$$ \varDelta I=\nabla \cdot \mathbf{G}, $$
(8)

solved in each iteration. The details of this method are available elsewhere, e.g. in [18].

4 Numerical results for considered cases

To examine the influence of the dendritic tip kinetics model on the rate of development of the undercooled region and columnar dendrites zone, the solidification of the binary alloy Al-4wt%Cu problem was solved in two geometries (Fig. 4). Configuration of boundary conditions ensures the 1D heat transfer and flat front (Fig. 4a) or 2D heat transfer and curved front (Fig. 4b). The initial temperature, T 0 is equal to 687 °C, and environment temperature, T env is 50 °C. Three cooling rates are considered, corresponding to heat transfer coefficients equal to 200, 500 and 1000 W/(m2 K).

Fig. 4
figure 4

Configurations of geometry and boundary conditions of considered cases

Full size image

Thermophysical properties of the considered alloy are listed in Table 1. Calculations were performed on a uniform mapped control volume (c.v.) grid. Numbers of divisions was 100 × 10 for the 1D heat transfer case and 50 × 50 for the 2D heat transfer case. The fully implicit Euler time integration scheme was utilized with time step equal to 0.005 s. At the end of each time step the Poisson equation for the indicator function was solved and equiaxed index was determined.

Table 1 Thermophysical properties of Al-4wt%Cu
Full size table

The developed computational model was verified with the built-in, simple model of solidification implemented in commercial software ANSYS Fluent. To obtain the full conformity of the both models, the linear relationship between volume fraction of solid phase and temperature was assumed. An excellent agreement was obtained. This computational model of binary alloy solidification was also positively verified for the Scheil solidification path model in [20, 21].

Temperature distribution in the considered domains, determined for the dendrite tip kinetics based on the dendritic tip stability parameter σ* = 0.02, and heat transfer coefficient 1000 W/(m2 K), is presented in Fig. 5. Both the solidus and liquidus isotherms are shown along with the columnar dendrite front position. They bound the regions where one of the two distinct grain morphologies is predominant. In the under-cooled liquid region, between liquidus isotherm (yellow line in Fig. 5) and the columnar dendrites front (green line in Fig. 5) equiaxed grains develop. In the considered model of equilibrium solidification, driven by the Scheil law, the solid phase appears at the liquidus temperature, so nucleation under-cooling is set equal to zero. The under-cooling liquid regions predicted in two considered cases are rather narrow, what corresponds to the simulation of Banaszek and Browne [19] performed for the same alloy and similar boundary conditions. A width of the under-cooled region increases in time while the temperature profile flattens so the conditions for equiaxed growth are more pronounced. Comparing the undercooled liquid zone at a chosen time (Fig. 5a and c) it is visible that it is more developed for the 2D case.

Fig. 5
figure 5

Temperature distribution in the two considered domains. Temperature maps (a) and (b) relates to 1D heat transfer 100 s and 200 s after the process start, (c) and (d) relates to 2D heat transfer after 100 s and 200 s. Position of the front is shown with a green line, liquidus and solidus isotherms are marked with thick yellow and black lines, respectively

Full size image

The evolution of the indicator function, predicted with the methodology given in [16, 19] (see section 3), is presented in Fig. 6 for the considered cooling rates and dendrite tip kinetics. The 1D and 2D cases are compared separately, for the selected heat transfer coefficient the impact of the dendrite tip kinetics on the equiaxed index is investigated. The equiaxed index is usually expressed as the function of a distance from the chill, e.g. [22]. This approach is convenient for cases where this distance is easily measured, namely in 1D geometries. For 2D and 3D cases the variation of the equiaxed index presented in [16] is adopted.

Fig. 6
figure 6

Equiaxed index as a function of time for: 1D heat transfer case (a, c, e) and 2D heat transfer case (b, d, f) for three different kinetics and various heat transfer coefficients: 200 W/(m2 K) (a and b), 500 W/(m2 K) (c and d) and 1000 W/(m2 K) (e and f)

Full size image

Regarding the 1D heat flow case the maximum in the indicator function is observed at the same time, roughly independent of the imposed kinetics. For the 2D heat transfer case the maximum is shifted towards the earlier times. In general, results for all numerical cases show that there is a large discrepancy between equiaxed predictors for FT numerical simulations using 2D and 3D dendritic growth models. Comparison between 1D and 2D results shows, the latter is more prone to CET, the maximum is shifted to earlier times.

Under-cooling averaged along the front (Fig. 7), after short time period stabilizes and approaches a virtually constant value, dependent on the used kinetics. The highest undercooling is observed for the dendrite tip kinetics based on the dendritic tip stability parameter σ* = 0.02. Undercooling is slightly more elevated for higher dimension of the domain. When the maximum temperature decreases below the liquidus temperature, the increase of undercooling at dendrite tips is observed.

Fig. 7
figure 7

Average under-cooling as a function of time for: (a) 1D heat transfer case and (b) 2D heat transfer case for three considered kinetics. Heat transfer coefficient is 500 W/(m2 K)

Full size image

A different criterion for the CET, called the constrained-to-unconstrained growth, was originally developed by Gandin [23] for the simplified model of 1D binary alloy solidification. The double front tracking approach was utilized, where two moving fronts across the domain representing the eutectic and columnar dendrite tips envelopes with prescribed growth kinetics were tracked. The mushy zone was bounded by isotherms of columnar dendrite tips and eutectic, so no nucleation and grain growth in the liquid undercooled zone appeared. The author concludes, that the maximum growth velocity appears at the CET position and the thermal gradient approaches zero in the in front of the columnar dendrite tips. The simulation outcomes (Fig. 8) show the temperature gradient does not achieve zero value but rather attains the roughly constant value. This is caused by the simple equilibrium model of solidification used in the analysis, so the solid phase appears on both sides of the front of columnar dendrite tips. Latent heat is released also in the undercooled liquid zone, so the local temperature raises and positive temperature gradient develops in front of the tracked envelope. This effect is more pronounced for 2D geometry than for the 1D one, and is observed for all analysed kinetics.

Fig. 8
figure 8

Averaged temperature gradient in front of the dendrite tip envelope for: (a) 1D geometry and (b) 2D geometry, for heat transfer coefficient 500 W/m2 K

Full size image

Additionally, the growth velocity of columnar dendrite tips (Fig. 9) related to the tip undercooling (Fig. 7) does not attain its local maximum, as was suggested by Gandin [23], but it rather tends to some constant value, dependent on the selected kinetics.

Fig. 9
figure 9

Average growth velocity of dendrite tips as a function of time for: (a) 1D heat transfer case and (b) 2D heat transfer case for three considered kinetics and the convective heat transfer coefficient equal to 500 W/(m2 K)

Full size image

5 Conclusions

The paper addresses the important issue of a proper choice of the stability constant in the commonly used dendrite tip kinetics models, which are an inherent part of computational simulations of alloy solidification, based on the cellular automaton, the front tracking approach or/and volume averaged single domain multiphase methods. The common practice is using the dendrite tip stability constant resulting from the marginal stability theory and experiments for 3D cases, despite of the computational model dimensionality. To analyze the correctness and accuracy of such an approach, 2D dendritic tip scaling parameter for aluminium-based alloy Al-4wt%Cu has been derived, then incorporated into the KGT dendritic growth model, compared with the original 3D KGT counterpart, and finally 2D and 3D versions of the modified Hunt’s analytical model for the columnar-to-equiaxed transition (CET) have been developed. This analysis shows significant differences in 2D and 3D versions of the stability constant, and the necessity to calculate sufficiently precise values of the tip scaling parameter in both cases.

The above is further confirmed through numerical calculations of the Al-4wt%Cu alloy solidification for three different dendrite tip kinetics by using the front tracking technique on a fixed control-volume grid. The results, including the mush-liquid front position, the calculated solutal under-cooling, the average temperature gradient at the dendrite tip envelope and the predictor of the relative tendency to form an equiaxed zone, are shown, compared and discussed.

In general, the analytical and numerical analyses show that there is a significant difference in the 2D and 3D representation of the CET effect, and the special attention should be paid for the proper selection of the stability parameter σ* - by taking into account the dimensionality of a dendrite tip growth model. In the case of 2D and 3D Hunt’s analytical map for the CET it is shown that, at the same the temperature gradient, the CET shifts to a higher velocity value of the 3D σ* stability parameter.