The boundary integral theory for slow and rapid curved solid/liquid interfaces propagating into binary systems

1. Introduction

The evolutionary behaviour of interfaces leading to developing multifarious pattern structures is of great significance in different phase transformation processes [1–5]. The Stefan-like moving boundary problems can be especially highlighted as very important tasks reflecting a complex dynamics of evolving patterns [6–12]. The phase interface in such problems is assumed to be atomically rough with thickness of the order of several interatomic distances. Not discussing this problem in detail (see works on phase-field modelling [13–26]) we shall neglect this atomic roughness and assume a sharp interface on the macroscopic length scale.

Methods of analytical and numerical description of the Stefan-type moving boundary problems can be classified as two problem-solving approaches: (i) to find the bulk temperature and solute concentration fields taking into account the phase interface conditions [9–11] and (ii) to formulate and solve the free-boundary problem for the phase interface function (deformation of the interfacial surface) [27–32]. The second approach can be directly used to evaluate the boundary integral suggested by Nash [27], discussed by Nash & Glicksman [28] and then independently detailed by Langer & Turski [29,30]. This method leans upon the Green function technique for the heat (and mass [33]) transfer equations, which allows a single integro-differential equation to be determined for the interface function. Note that Langer & Turski [29,30] developed their original theory to obtain the boundary integral for an Ising-like symmetric model when two phases contact each other at the interface. This theory has become the basis for stability and selection problems of dendritic interfaces [34–45]. Until quite recently, such a coupled thermo-diffusion problem with the help of boundary integral theory was studied only within the framework of the linear temperature gradient model by Kessler et al. [46]. Generalization to the complete coupling of the fields on the interface has been treated in [33].

The boundary integral method has a great impact on interface evolution and pattern formation in different free boundary problems (e.g. with the selection mode of stable dendritic growth [37–39]). This method can also be applied for modelling and testing several crystal growth theories or for comparison of the model results with the phase-field simulations. We demonstrate below how to use the boundary integral theory to also calculate interfacial shapes. Concluding this introductory section, let us especially emphasize that the present work generalizes the boundary integral theory for parabolic heat and mass transfer problems, extends it to the hyperbolic mass transfer in liquid, indicates how to simplify it in the case of steady-state interface propagation and demonstrates how to obtain the classical Ivantsov and Horvay–Cahn stationary solutions and growing shapes. Patterns selected in quasi-stationary solidification are modelled by the boundary integral and numerically solved by a specially developed iteration technique.

2. Thermo-chemical interface propagation under local equilibrium conditions (parabolic heat and mass transport)

Consider the solid/liquid interface propagation into an undercooled binary mixture caused by temperature and solutal gradients when the anisotropies of surface energy and growth kinetics play a preferable role in the evolutionary behaviour of the interface function (deformation of the interfacial surface). To obtain a time-dependent integro-differential equation for this function let us use the Green function technique previously developed by Nasha & Glicksman [27,28] and Langer & Turski [29,30]. A solution procedure is divided into two simple steps which are concerned with a one-by-one consideration of purely thermal and thermo-chemical growth regimes.

(a) Thermally controlled interface propagation

First of all, we consider the thermally induced motion of the solid/liquid interface into an undercooled pure liquid. Let initially the solid/liquid interface be flat in the x-plane, which is perpendicular to the spatial direction z. The interface becomes curvilinear in the course of time due to the presence of instabilities and preferable local growth directions (figure 1). This curvilinear interface is designated by means of the interface function z_interface=ζ(x,t), where t is the phase transition time.

• Download figure
• Open in new tab
• Download PowerPoint

To find the interface function equation, two possible cases can be applied: (i) equal thermal diffusivities D_T in both phases [37,38] and (ii) a one-sided phase transition model where the transport flux is negligible from one of the phases at the moving interface [33]. Let us now consider the first of these approximations. The second one is used below when analysing the chemically controlled interface growth.

At the moving phase transition interface ζ(x,t), the Gibbs–Thomson and heat balance conditions hold true,

and where T(x,z,t) is the temperature, is the anisotropic coefficient describing the kinetic growth mode, T_f is the phase transition (freezing) temperature of the planar front at , d_c is the anisotropic capillary length, Q is the latent heat of phase transition released per unit volume of the solid, c_p is the specific heat, V is the constant growth velocity, ds is the vector element of the surface area pointing to the liquid phase, subscripts ‘s’ and ‘l’ designate the solid and liquid phases, respectively, and the interface curvature K can be represented in the two- and three-dimensional cases in the form

Introducing the coordinate system moving with velocity V , one can write down the heat transport equation in the liquid and solid phases as

Green’s function G(p|p₁) corresponding to equation (2.4) satisfies the following equation:

where δ is the Dirac delta function, p is a short designation for the space–time point (x,z,t), the operators on the left-hand side act on the subscripted variables x₁, z₁, t₁ and G vanishes for t₁>t.

Multiplying equations (2.4) and (2.5) by G and T, respectively, with their further subtracting leads to

where Λ₁ is a region that includes the point p and excludes the interface, S₁ is the overall bounding surface of Λ₁ and ds₁ is the vector element of the area on S₁ pointing along the inward-directed normal to this surface.

Green’s function determined by equation (2.5) can be expressed in terms of the corresponding Fourier transformation as

where r is the vector (x,z), the contour C goes from to and passes over all singularities in the ω-plane in order to satisfy the causality condition [30], i is the imaginary unit, and parameter j equals 2 and 3 in the two- and three-dimensional cases, respectively. The integral kernel entering equation (2.7) can be represented as The inverse Fourier transform of equation (2.8) gives

Because the temperature field is continuous at the phase transition interface, the contribution proportional to ∇₁G in equation (2.6) cancels out. Now, we are able to calculate the second integral term containing ∇₁T in equation (2.6) on the basis of expression (2.2). Keeping in mind that the last term in equation (2.6) is equal to [30], we get

Here, stands for the temperature far from the solid/liquid interface and ζ₁ designates ζ(x₁,t₁).

If we now move the point p closer to the phase interface and take into account equations (2.1) and (2.10), we come to

where is the dimensionless undercooling and is the kinetic coefficient, and

Let us introduce the following dimensionless variables (ρ is a characteristic length scale):

Substituting (2.12) into the integral equation (2.11) and omitting primes for simplicity, we arrive at where in the two-dimensional case and in the three-dimensional case. Here, P_T=ρV/(2D_T) is the thermal Péclet number, and K(ζ(x,t)) is defined by expressions (2.3). Equation (2.13) represents the dimensionless integro-differential equation for the interface function in two and three spatial dimensions if the phase transition is thermally controlled.

An important point is that expressions (2.13) and (2.14) represent the basis for the previously developed stability and selection theories of dendritic growth [39–43] in the case of two-dimensional geometry. In addition, equations (2.13) and (2.15) transform to the corresponding expressions of [2,37] in the limiting case ∂ζ/∂t=0, and in the absence of growth anisotropy if we are dealing with the three-dimensional case. Note that equations (2.13)–(2.15) have been used by Barbieri & Langer [38] in the absence of interface atomic kinetics, in two- and three-dimensional growth geometries.

(b) Thermally and chemically controlled interface propagation

In this section, we consider the so-called thermo-chemical interface propagation occurring as a result of temperature and solutal gradients when the anisotropies of surface energy and growth kinetics dictate the preferable growth direction. Such a dynamic scenario of interface propagation occurs in a binary mixture during the intermediate stage of phase transformation between the diffusion-limited and pure thermally controlled regimes [47,48].

An isobaric binary system consisting of non-interacting A- and B-kinds of particles (solvent and solute) is under study when the phase transition occurs without the mass density change (i.e. without shrinkage of phases). We also neglect any convective processes of either heat or solute (A- and B-particles) and, as above, consider the conductive mechanism of heat (and mass) transfer. Keeping in mind that the chemical diffusion in the solid phase is much slower than that in the liquid phase [9,49–53], we consider the diffusion of B-particles in liquid only (one-sided model). Therefore, the mass balance equation, written out in addition to the heat balance equation (2.2), takes the form

where D_C is the chemical diffusivity of the solute (B-particles), k₀ is the equilibrium partition coefficient, C_l is the solute concentration of the binary mixture in the liquid phase and C_i is the interfacial concentration. The Gibbs–Thomson condition (2.1) for thermo-chemical growth has the form where m designates the liquidus slope. The mass transport equation in the moving frame of reference is governed by

It is important to note that the interface kinetics previously introduced in the thermal and thermo-chemical problem of solidification has been validated independently by means of other alternative approaches [54,55]. This makes possible a flexible description of thermal and chemical processes at different stages of dendrite growth [56].

Equations (2.4) and (2.18) are equivalent and coincide after replacing D_T by D_C. Therefore, by analogy with the solution (2.6) for temperature, the solution for the mass transport equation is described by

Here, the integral of ∇₁[C_l(p₁)G_C(p|p₁)] becomes zero in accordance with the Kelvin–Stokes theorem of integration of the nabla operator about a closed path S₁. Next, one can easily find the integral corresponding to the contribution G_C(p|p₁)∇₁C_l(p₁) by means of the boundary condition (2.16) if neglecting the concentration flux in the solid phase. The last contribution in (2.19) equals the solute concentration in liquid far from the solid/liquid interface (for details, see [30]). As a result, we are able to rewrite expression (2.19) as Combining expressions (2.13) and (2.17), one can obtain the interfacial concentration

Now the interface function ζ(x,t) is determined from expressions (2.20) and (2.21) and satisfies the following time-dependent integro-differential equation:

where in the two-dimensional geometry and in the three-dimensional geometry. Here, P_C=ρV/(2D_C)=P_TD_T/D_C stands for the chemical Péclet number. Note that P_C in equations (2.23) and (2.24) are pulled out of the integrals because they include only constant values of V and ρ related to the quasi-stationary growth mode.

In consideration of the presently formulated thermo-chemical boundary integral, we have used a simple model of a dilute binary mixture. A more realistic boundary integral model of multi-component non-isothermal mixtures can be developed if the diffusion fluxes are written via chemical potentials and, respectively, thermodynamic potentials from all elemental and energetic contributions of the whole mixture, as has been done in models based on the phase-field method [57,58].

Concluding this section, let us underline that the boundary integral that describes the curvilinear interface function ζ(x,t) for the thermo-chemical Stefan problem is given by equation (2.22), whose thermal and concentration contributions are determined by expressions (2.14), (2.15), (2.23) and (2.24). As a special note, the isothermal solute diffusion-limited solidification is defined by equation (2.22). Namely, setting and P_T→0 (in this limiting case, determined by equations (2.14) and (2.15) vanishes), this case is also described.

3. Thermo-chemical interface propagation under local non-equilibrium conditions (parabolic heat and hyperbolic mass transport)

The aforementioned thermo-chemical boundary integral theory can only be applied for the description of local equilibrium phase transition. Modern experimental facilities allow us to obtain large undercoolings and cooling rates leading to substantial temperature and concentration gradients as well as high rates of phase interface propagation. For instance, the electromagnetic levitation technique allows the liquid phase to be undercooled by up to 400 K below its liquidus temperature. This leads to the highest solidification velocities ranging from 10 to 100 m s⁻¹ [59]. A theoretical description of such deeply undercooled phase transitions depends upon the hyperbolic mass transport equation for solute diffusion [47]. Such a hyperbolic approach has been validated by comparison with experimental data [60], by atomistic simulations of the solute trapping effect for fast fronts [61] and by coarse-graining derivations of the phase-field equations [62].

(a) Mathematical background of the hyperbolic mass transport

Taking the aforesaid into account, we consider the phase interface ζ(x,t) propagating into a metastable medium under local non-equilibrium conditions which are described by the hyperbolic-type mass transfer equation [47] in the laboratory coordinate system

where τ_D stands for the local diffusion relaxation time of mass flux to its steady-state value. As the temperature field relaxes to its local equilibrium much faster, we use, as before, the parabolic temperature conductivity equation (2.4).

At the curved solid/liquid front ζ(x,t), the mass conservation has the form

where the velocity-dependent chemical segregation coefficient is [63] where v=V +∂ζ/∂t, and V_D and V_DI are the bulk and interface diffusion speeds, respectively. Note that the thermal fluxes at the phase interface, as before, are connected by means of the boundary condition (2.2).

The Gibbs–Thomson condition now becomes

where m_v=dT/dC_l is the change of temperature at the interface with respect to the change of concentration as a function of the interface velocity at constant pressure and T_Q=Q/c_p. The velocity-dependent liquidus slope m_v takes the form [64]

(b) Boundary integral equations

The Green function G_CH for the hyperbolic mass transport equation (3.1) satisfies

where subscript CH designates the hyperbolic solution for concentration. For two and three dimensions, G_CH has been determined by Morse & Feshbach [65] as where and where , j=2 and j=3 in the two- and three-dimensional cases, respectively, ρ is a characteristic length scale, J₁ is the Bessel function of the first kind, θ is the Heaviside step function, τ_*=2τ_DV/ρ, , τ=t−t₁, W=V/V_D, υ=2D_C/V_D, P_C=ρV/(2D_C)=W²/τ_* is the solutal Peclét number, and all length and time scales, as before, are measured in units of ρ and ρ/V , respectively. Note that Green’s functions (3.8) and (3.9) turn to zero at τ<R_2D and τ<R_3D, respectively. An important point is that solution (3.7) transforms to the corresponding parabolic Green’s functions determined in the preceding section in the limiting case τ_*→0 (see also [27,33]).

We especially note that Morse & Feshbach [65] have not considered transport velocities reaching and exceeding the characteristic speeds for the propagation of disturbances, i.e. they have written the Green functions (3.8) and (3.9) only for the case of velocities W=V/V_D<1. In the case of alloys, the solidification front can reach or even exceed the atomic diffusion speed. Therefore, we have included in equations (3.8) and (3.9) the conditions R_jD<τ and R_jD≥τ, which are equivalent to regimes with the growth velocity smaller than the maximum diffusion speed, W<1; V <V_D, and regimes with the growth velocity larger than the maximum diffusion speed, W≥1; V ≥V_D, respectively.

Now, the interfacial concentration at z=ζ should be written as

and where , k_v is calculated at the time moment t−τ and x₁ must be replaced by x₁ in the two-dimensional case. Note that P_C in equation (3.11) are pulled out of the integrals because they include only constant values of V and ρ related to the quasi-stationary growth mode. Two inequalities in equation (3.12) follow from equations (3.8) and (3.9) as well as from the conditions of the diffusionless solidification C_i=0 and k_v=1 with V ≥V_D (W≥1; see equation (3.3)).

A unified integro-differential equation for the interface function ζ(x,t) is determined from the Gibbs–Thomson condition (2.17) or (2.22) with allowance for the hyperbolic solutions (3.10)–(3.12) as

This equation describes the interface evolution for the thermo-chemical Stefan model when the thermal field is defined by the solution of the differential equation of the parabolic type and the solute field is described by the differential equation of the hyperbolic type. As follows from equations (3.10)–(3.12), the complete equation (3.13) is true for V <V_D. If V ≥V_D, one obtains . As a limiting case, equation (3.13) transforms to equation (2.22) for the parabolic thermo-chemical model when V ≪V_D, W≪1.

4. Steady-state interface propagation

The boundary integrals derived in the preceding sections can be simplified in the practically important solidification process with a constant velocity V , ∂ζ/∂t=0. This section is concerned with modifications of the thermal and concentration integrals for parabolic and hyperbolic transport.

(a) Thermo-solutal parabolic problem

The right-hand side of equation (2.14) can be evaluated over τ for two dimensions, ζ=ζ(x), by taking into account that the modified Bessel function of zero order is expressed as [66]

Then, equations (2.13) and (2.14) have the form of a single modified boundary integral, which determines the curved interface function ζ(x) for the two-dimensional thermal problem.

Keeping in mind that the modified Bessel function of -order [66] reads as

integrating the right-hand side of expression (2.15) over τ and substituting the result into equation (2.13), one can obtain the steady-state boundary integral in the three-dimensional case, ζ=ζ(x,y), as This integral characterizes the interface function ζ(x,y) for the thermal problem in three dimensions. Note that equations (4.1) and (4.2) have been previously used by Pelcé [2] and Kessler et al. [46] in their analyses of stable dendritic growth.

In the case of the coupled thermo-solutal problem, keeping in mind equations (2.21) and (2.22), we have the following equation for the interface function, which includes the thermal and concentration contributions and :

Here, is determined by equations (4.1) and (4.2) in the two- and three-dimensional case, respectively.

Taking into account a certain analogy between the thermal and concentration problems, we are now able to find the concentration integrals as

in the two-dimensional space and in the three-dimensional space.

Concluding this subsection, let us especially underline that equation (4.3) determines the steady-state interface function ζ(x,y) for the thermo-solutal parabolic Stefan problem, and the integrals and are determined by expressions (4.1), (4.2) and (4.4), (4.5).

(b) Thermo-solutal problem in the case of hyperbolic mass transfer

Let us now integrate expressions (3.11) for the quasi-stationary solidification. Completing the squares in the arguments of the hyperbolic cosine and Bessel function J₁ in expressions (3.8) and (3.9) and integrating over τ, we arrive at the modified concentration integrals in the two- and three-dimensional geometries, respectively [67]:

and Here, K₀ and K_1/2 represent the modified Bessel functions, and b_H=(1−P_Cτ_*)|x−x₁|²+(ζ(x)−ζ(x₁))².

Expressions (4.6) and (4.7) are true for V <V_D (as in the non-steady case described by equations (3.10)–(3.12)). If V ≥V_D, one obtains in equations (4.6) and (4.7). This result has a clear physical meaning [68]: when the interface velocity V is equal to or greater than the diffusive speed V_D, a profile of concentration cannot be created ahead of the interface and the system solidifies by the diffusionless mechanism. Note that the transition to diffusionless solidification proceeds sharply with an abrupt drop point at V =V_D in the kinetic curves ‘dendrite velocity–undercooling’ [47].

The steady-state parabolic heat transfer therewith is described by integrals from expressions (4.1) and (4.2), whereas the interface function ζ(x,y) follows from equation (3.13) at ∂ζ/∂t=0. As a result, the equation for ζ is described by

(c) Selection criterion of stable dendritic growth

This subsection is concerned with a selection criterion σ*=2d₀D_T/(ρ²V) that determines a combination between the dendrite tip diameter ρ and its tip velocity V . We obtain this combination on the basis of dendritic growth theory developed in a similar paper in the present theme issue [20]. The criterion for σ* given below generalizes previously known criteria and comprises the hyperbolic mass transfer mechanism in a binary melt, growth kinetics and forced convection. The main idea of how to derive this criterion consists in developing the laws for formal transition between the hyperbolic and parabolic equations that describe the solute concentration in the liquid phase.

Let us initially emphasize an important feature of hyperbolic integrals (4.6) and (4.7) characterizing the steady-state dendritic growth at high solidification velocities V (at high growth Péclet numbers P_T). This feature resides in the fact that the integrals (4.6) and (4.7) (for the transport described by the hyperbolic equation) transform into equations (4.4) and (4.5) (which are true for the transport described by the parabolic equation) after introducing new spatial coordinates (subscript ‘new’) of the form

(and , in the three-dimensional geometry). Indeed, keeping in mind that (and in the three-dimensional case), and k_v=k_v becomes k₀, one can conclude that expressions (4.6) and (4.7) become identical to expressions (4.4) and (4.5). Expressions (4.9), in turn, imply that the hyperbolic problem transforms into the parabolic one after substitution by ρ_new in coordinates x_new, y_new and z_new. It is important to highlight here that the selection parameter σ* (or combination ρ²V) remains constant during dendritic growth. By this we mean that the tip velocity (1−W²)V should be replaced by V_new (ρ²_newV_new=const.), while the chemical Péclet number P_C=ρV/(2D_C) becomes .

Expression (3.25) of [20] determines the selection criterion for the parabolic mass transfer mechanism. Taking the aforementioned remarks into account, we rewrite this criterion for the hyperbolic mass transfer mechanism and n-fold symmetry of crystal growth as (for the sake of simplicity, we omit subscript ‘new’ and designate through V and ρ the dendrite tip velocity and its tip diameter for the hyperbolic problem)

where d₀ is the capillary constant, U is the velocity of forced flow, α_d is the surface energy anisotropy, β₀ is the kinetic constant, and σ₀, b, a₁, a₂ and δ₀ are constants.

Here, k_v=k_v(V), m_v=m_v(V) and C_i are determined by relations (3.3), (3.5) and (3.10), respectively. Let us note that, in criterion (4.10), we have rescaled only the chemical Péclet number (the growth Péclet number P_g remains unchanged). This is due to the fact that the rescaling procedure is concerned only with the chemical dendrite selection criterion (the second term in curly brackets of (4.10)) induced by changing the mass transfer from parabolic to hyperbolic. The selection parameter (4.10) has only the thermal contribution in the rapid solidification limit W≥1 (V ≥V_D) when the concentration distribution in liquid is frozen and does not contain any inhomogeneities. The hydrodynamic contribution entering in completely coincides with the parabolic growth problem because this contribution plays a key role only in the limit of slow velocity V when the hyperbolic mass transfer does not occur and the crystal growth proceeds in accordance with the parabolic heat and mass transfers.

The criterion (4.10) contains previously found parabolic selection criteria discussed in [20,69]. The criterion describes the thermally and solutally limited regimes under the control of the anisotropy of surface energy. Limiting cases of the solely kinetic regime of dendritic growth, an influence of kinetic anisotropy and the transition from the thermo-solutally controlled to the kinetically limited regime are given in [43,56]. In these works, asymptotic solutions and sewing with the kinetic regime are shown and discussed. Using asymptotics for the kinetic regime, one can obtain the selection criterion, which includes the kinetic anisotropy additionally to equation (4.10).

In summary, the criterion (4.10) extends the growth criteria derived for the transport described by the parabolic equation and includes hyperbolic mass transfer playing a substantial role in the case of rapid dendritic growth processes [47]. This criterion also generalizes the described limiting regimes and gives a stable relation between V and ρ for the n-fold crystal growth symmetry with allowance for the kinetics and surface energy anisotropies and a forced convective flow.

5. Steady-state Ivantsov solutions

This section is concerned with the question of how to derive the steady-state solutions describing the growth of a parabolic platelet (two-dimensional growth) and a paraboloid of revolution (three-dimensional growth) in the case of parabolic and hyperbolic transfer. At the first instance, we turn our attention to the classical Ivantsov solutions for the thermo-chemical problem.

(a) Ivantsov solutions

Let us demonstrate that the boundary integrals (2.13)–(2.15) and (2.22)–(2.24) include the classical Ivantsov solutions [70,71], which describe the steady-state propagation of a parabolic platelet and a paraboloid of revolution. We assume that the surface energy and kinetic effects are negligible (d_c=0, β=0) and ∂ζ/∂t=0.

Let us begin with the thermally controlled growth of a parabolic platelet and consider expressions (2.13) and (2.14). Initially, we substitute the parabolic shape ζ(x)=−x²/2 and replace the integration variable τ by ω as

Then we integrate the right-hand side of expression (2.14) over x₁ and conclude that the resultant expression is independent of x. Combining the obtained expression and equation (2.13), we come to the liquid undercooling Δ, which coincides with the thermal Ivantsov solution [33,40,46,70,71] Here, we use the following integral expression [67]:

When we are dealing with the three-dimensional motion of a paraboloid of revolution ζ(x,y)=−x²/2−y²/2, we replace two variables τ and y₁ by ω and z₁, respectively, in equation (2.15) as

Integration over x₁, and then over z₁, gives Here, we take into account that [67]

The next step is to differentiate and integrate with respect to its parameter with allowance for

As a result, we arrive at the Ivantsov solution describing the liquid undercooling Δ for a paraboloid of revolution [33,70,71],

Combining the obtained equations (5.1) and (5.5), we are able to write down the generalized solution that is valid for the two- and three-dimensional cases

where for two-dimensional and j=1 for three-dimensional geometries.

In this case of the thermo-chemical Ivantsov solutions, we formally replace P_C by P_T in integrals (2.23) and (2.24) and keep in mind the concentration factor (1−k₀)C_i. Hence, we get

By means of further substitution of into equations (2.21) and (2.22) at d_c=0, β=0 and ∂ζ/∂t=0, we arrive at the familiar interfacial concentration [71]

Expressions (5.1)–(5.8) describe the thermo-chemical propagation of a parabolic platelet or paraboloid of revolution with a constant velocity (for details, see [72]). The aforementioned equations for interface motion (2.13)–(2.15) and (2.22)–(2.24), as a consequence, contain the steady-state Ivantsov solutions so that

Finally, the temperature and solute diffusion fields in liquid around the growing parabolic platelet (two dimensions) and paraboloid of revolution (three dimensions) are

(b) Solutions for the hyperbolic equation of mass transport

Let us now consider the question of how to obtain the steady-state interfacial concentration at the parabolic platelet and paraboloid of revolution in the case of hyperbolic mass transport. This problem was considered in [68,73]. This derivation leans upon equations (3.10), (4.6) and (4.7) and the following auxiliary expression [67]:

Again, in the two-dimensional case we introduce the new variable ω=(x−x₁)/α, integrate over x₁ and arrive at the integral, which is independent of x. In the three-dimensional case, x=(x,y) and x₁=(x₁,y₁); we additionally introduce the second variable z₁=(y−y₁)/(x−x₁), and integrate over x₁ and z₁. Omitting mathematical details, we have

where and j=1 for the two-dimensional and three-dimensional geometries, respectively.

The interfacial solute concentration can be found from expressions (3.10) and (5.11) as [47,68]

The solid/liquid interface therewith moves in accordance with expression (4.8), where the integrals and are determined in (5.6) and (5.11).

The concentration integrals (4.6) and (4.7) at P_Cτ_*=W²=1 become zero and the Green functions (3.8) and (3.9) turn to zero in the case W>1. This implies that and in accordance with expression (3.10). By this we mean that the concentration profile is constant (diffusionless solidification) ahead of the phase interface when the interface velocity V is equal to or greater than the diffusive speed V_D [47]. Expressions (5.12) coincide with those previously found by means of the Ivantsov and Horvay–Cahn methods [68,73] for hyperbolic mass transport. Solutions (5.11) and (5.12) are correct only in some area behind the dendrite tip. Far from the tip, in the region where the undercooling may be much smaller than at the tip, the dendrite velocity can be lower than the diffusion speed, V ≪V_D, W≪1. In such a case, the dendrite may grow under local equilibrium conditions and solutions (5.11) and (5.12) transform to the classical Ivantsov solutions obtained using the parabolic transport equations [70,71].

6. Steady-state Horvay–Cahn solutions

A number of experimental data show that in many cases the growing shapes of dendrites are non-symmetric. For instance, ice crystals evolving in pure water and in water–ethyleneglycol solution or dendritic evolution in aluminium nitride, lithium niobate or lithium tantalate represent several examples of the sixfold symmetry which is caused by the hexagonal crystalline anisotropy [4,5,74–78]. In these cases, dendritic shapes do not satisfy the Ivantsov solutions [70,71] and one must use the Horvay–Cahn solutions characterizing possible lower symmetric dendritic shapes [79–81]. In this section, we demonstrate how to derive their analytical solutions from the boundary integral theory.

(a) Parabolic heat and mass transfer equations

Horvay & Cahn [79] showed that the shape of a growing dendrite can be an elliptical cross section, such that the dendritic surface in dimensionless variables is given by

where 2D_T/V =ρ/P_T is the length scale and ρ represents the dendrite tip diameter of any of its cross sections (or the average value of the dendritic tip curvature radii). The dendrite interface is described by ω=P_T and a represents the aspect ratio of the elliptical cross section of the dendrite. The cross section is circular if a=0 (|a|<P_T). Equation (6.1) represents a symmetrized case of the corresponding Horvay–Cahn equation, keeping in mind the dimensionless radii of curvature P_T−a in the x–z plane and P_T+a in the y–z plane. Thus, the average dimensionless radius of curvature is P_T [80].

In dimensional form, expression (6.1) becomes

where ω_d=ωρ/P_T and a_d=aρ/P_T, and the subscript d designates the dimensional variables and parameters. The dimensionless interface function and undercooling are determined by expressions (2.13) and (2.15).

In the steady-state propagation of an elliptical paraboloid into an undercooled liquid with constant velocity V (∂ζ/∂t=0), equation (2.15) reads

This equation can be rewritten by means of expression (6.2) as

Now, replacing the integration variables τ and y₁ by ω₁ and z₁, respectively, using equation (5.2), we arrive at Integrating the right-hand side of expression (6.5) over x₁, we get The next step is to integrate the right-hand side of equation (6.6) over z₁ taking into account the integral expression (5.4) as

Replacing variable ω₁ by , we have

where

Using the method of differentiation and integration of J(ω_d), we have

At the interface, ω_d=ρ, we get where and Expression (6.13) for the melt undercooling gives a relation between the unknown dendritic parameters ρ and V . The present solution (6.11) transforms to the Horvay–Cahn twofold solution [80] after transition to the dimensionless parameter a=a_dP_T/ρIt is important to note that the case of the circular cross section, a=0, gives the Ivantsov solution [70,71].

The dimensional Horvay–Cahn temperature distribution in liquid has the form [79–81]

The Horvay–Cahn-type solution for the thermo-chemical problem can be found by analogy with the temperature problem (6.3)–(6.15). Using similarity between the thermal and concentration problems, we have

and

The function ω=ω(x,y,z) is given by equation (6.1). This equation has three roots, but only one of them is positive for all possible values of parameters. This root is found from the Cardano formula

where and

Analysing the Horvay–Cahn solutions [80], it is suitable to use the parabolic coordinates (ξ,η,φ) connected with the Cartesian ones (x,y,z) as

The dendritic shape steadily growing in the z-direction as an elliptical paraboloid is described in the parabolic frame of reference by This expression determines ω as a function of (ξ,η,φ) in the form where and

Thus, the Horvay–Cahn solutions completely follow from the boundary integral theory. The thermo-solutal parabolic problem within the Horvay–Cahn formalism is described by expressions (2.21) and (6.14)–(6.17).

(b) Parabolic heat and hyperbolic mass transfer equations

Replacement of the variable z by z_h,

in the parabolic solute transport equation leads to the hyperbolic equation taking into account the finite speed V_D of atomic diffusion [73]. Defining , we get the average tip diameter, chemical Péclet number and elliptical paraboloid equation (6.2) as

The integral equation of the interface motion in dimensionless form is given by equation (4.8), where reads as (see equations (4.6) and (4.7) and the discussion around them)

Further, in this section we consider only the case of velocity which is smaller than the diffusion speed, i.e. W<1, V <V_D.

Returning to the dimensional coordinates and assuming that the elliptical paraboloid does not change its shape and moves with constant velocity V (i.e. ), we get from equation (6.28) that

Replacing two integration variables τ and y₁ by ω₁ and z₁, respectively, defined by equation (5.2), we obtain Now, replacing the integration variable x₁ by and integrating the right-hand side of equation (6.30) over , we get Integrating the right-hand side of equation (6.32) over z₁ and taking into account equation (5.4), we evaluate the chemical integral (6.31) as Replacing ω₁ by , we obtain where Further, applying the method of differentiation and integration of J_CH(ω_d), we get At the interface, ω_d=ρ, one obtains where

Expressions (6.37) and (6.38) transform to equation (5.11) in the limiting case a_d=0 of the circular parabolid. Note that integral (6.37) is identical to integral (6.16) in the case of the parabolic mass transfer mechanism when V ≪V_D, P*_C→P_C and a tends to a_dP_C/ρ. The concentration distribution in the present case is given by expressions (6.17) when P_C is replaced by P*_C, whereas the temperature field is governed by equations (6.15).

As a final note of this section, the boundary integral theory includes the Horvay–Cahn steady-state temperature and solute concentration solutions within the framework of parabolic and hyperbolic transport equations.

7. A globular shape of dendritic tips at high Péclet numbers

In this section, we demonstrate that the solid/liquid interface, which is determined by the boundary integrals (2.13), (2.22) and (3.13), has a circular (spherical) shape in the steady-state solidification conditions with high Péclet numbers in the case of either the parabolic or hyperbolic description of transport. Setting and , we conclude that integrals (2.14), (2.15), (2.23), (2.24) and (3.11) tend to zero, whereas conditions (2.13), (2.22) and (3.13) for the interface function become

where the interface curvature K is given by expressions (2.3).

Substituting K for the two-dimensional case into (7.1), we get

This equation can be easily integrated in the form where ζ₀ and C₀ are arbitrary constants. Equation (7.2) shows that the interface function ζ(x) represents a branch of the circle. This can be easily demonstrated by means of the squaring operation and transformation to the canonical form.

In the three-dimensional case, we have from equations (2.3) and (7.1)

where i and j are the unit coordinate vectors, a_x, a_y, C_x, C_y are constants and a_x+a_y=−B. This equation implies Now, integrating equation (7.3) and taking into account expression (7.4), we obtain the surface with a_x=a_y=−B/2 as Equation (7.5) describes a spherical surface similar to that of the two-dimensional geometry (7.2). Asymptotically, the boundary integral equations define the growth of spherical surfaces (or nuclei), which is caused by the constant driving force entering into the right-hand side of equation (7.1).

Now, we solve equations (7.1) and (7.5) to check the shape of the dendrite tip in the anisotropic and isotropic cases. Because the anisotropic capillary length d_c and anisotropic kinetic growth coefficient β in the three-dimensional case have the following form (see [48] and references therein):

and it is convenient to transform Cartesian coordinates to spherical coordinates defined by where α_β is the kinetic anisotropy parameter.

To obtain the growth velocity V and the tip diameter ρ in equations (7.1) and (7.5), we use the selection criterion σ* as

Because equation (7.1) has been found within the limits and and we have to obtain the surfaces with the finite and non-zero values of V and ρ from equation (7.9), we use σ* from the kinetic growth mode [56], for which we merely take the large and finite value of the Péclet number P_T=1. In equation (7.10), we have and . The constants a₁, a₂ and δ₀ are given by

For the computation of shapes by equations (7.1) and (7.5) using equations (7.9)–(7.14), we have the chosen material parameters of an Ti-Al alloy given in table 1. The calculations have been made for the steady-state growth at the dimensionless undercooling Δ=1. The results of calculations are shown in figure 2, which demonstrates the areas around the dendrite tips in the anisotropic (a) and isotropic (b) cases. At the large value of the Péclet number, P_T=1, chosen in calculations, these shapes are qualitatively consistent with the globular shape appearing at the high-growth regimes of crystal formation [83–87]. Computations by the phase-field model of Karma & Rappel [88] showed that the crystal growth morphology also exhibits non-branching shapes of a globular form at very high undercooling [89]. Even though the characteristic scales of patterns of figures 2 and 3 are different (d₀≪10D_T/V), these crystal shapes have a clear qualitative similarity: the deformed semi-sphere (figure 2a) has a similarity to the anisotropic spherical globe (figure 3a), and the spherical tip of the dendrite (figure 2b) is similar to the isotropic spherical globe (figure 3b). From this it follows that the solution (7.5) describes the globular tips of dendrites qualitatively consistent with the results of phase-field modelling for high growth rates of crystals.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

With regard to the comparison of patterns reproduced by the boundary integral method (figure 2) and the phase-field model (figure 3), a natural question arises about the preferable applicability of these approaches depending on the modelled system and length scales. Indeed, the phase-field model is a well-established thermodynamically consistent method, which is nowadays developed for many multi-component and multi-phase systems (see [57] and references therein). It must be applied to the transformations and reactions in which the diffuse interface between different states plays a specific role [90,91]. In solidification processes, the diffuse interface has a width of several lattice parameters and a whole pattern, which can be computationally reproduced by the phase-field model, and has a scale from several microns to fractions of millimetres. The boundary integral method is based on the sharp interface formulation in which the thickness of the interface is negligibly small, i.e. it has zero width. Owing to this feature, the boundary integral method may cover macro-scales (up to metres in lengths), for instance in modelling of dendrites and grains growing or melting in the Arctic sea. In addition to this, it can reproduce the dendritic pattern much faster than the phase-field method. Therefore, patterns modelled by the boundary integral method may be used as the first benchmarks for the phase-field model, which can then be used for more precise solutions.

Concluding this section, one has to remark that equation (7.1) represents the Gibbs–Thomson condition in the limit of infinite Péclet numbers. A consistency of such a limiting case with the velocity and acceleration-dependent Gibbs–Thomson equation [92] would be a subject of further analysis.

8. Computational modelling

(a) Model equations

Neglecting the concentration gradient in solid, ∇C_S=0, let us introduce concentration , where C_i(r) is the concentration at the interface and is the concentration of a solute at infinity. All analytical treatments as well as all numerical approximations are made in the dimensionless coordinates

where h_D=D_C/V_D is the length scale for inter-diffusion.

We shall model a two-dimensional dendrite shape under steady conditions with constant interface velocity. Then, the diffusion mass transport in the region Ω (figure 4) in coordinates (8.1) is described by

Substitution of the concentration (3.10) with the integral (4.6) into the mass balance (3.2) gives the equation for the steadily propagating interface, where the normal n is directed from the solid to the liquid. Then, the boundary integral for the steady interface takes the following form: with ζ_i the whole interfacial contour, including the phase interface ζ_SL and boundaries of the computational domain shown in figure 4. In equation (8.4), we have the Green function of the operator (because one gets L⁺g(r−r′)=−δ(r−r′)), and the kernel. Here, K₀ and K₁ are the modified Bessel functions of the zero and first order, respectively, is the norm in the two-dimensional domain, and is the modified Péclet number for the inter-diffusion, which depends on the ratio W=V/V_D between the interface velocity V and the maximum diffusion speed V_D.

• Download figure
• Open in new tab
• Download PowerPoint

The boundary integral (8.4) is obtained for W=V/V_D=const, and it is considered as the equation for the normal gradient of concentration, ∂u(r)/∂n, which drives the interface. Equation (8.4) and expressions (8.5)–(8.8) are true for an interface velocity smaller than the diffusion speed, i.e. at W<1, V <V_D. If the interface moves with a velocity equal to or greater than the solute diffusion speed, W≥1, V ≥V_D, diffusionless growth with , u(r)=0 occurs [60,63], and both sides of the boundary integral equation (8.4) are reduced to zero (see also steady-state solutions (4.6) and (4.7) and the discussion around them). A boundary integral that is consistent with the analytical results does not exist [68]: for the case W≥1, V ≥V_D, an arbitrary crystal shape is possible in the isothermally solidified binary mixture. Finally, in the local equilibrium limit, W≪1, V ≪V_D, equations (8.4)–(8.8) are reduced to the previously known boundary integral of Saito et al. [31], obtained for the analysis of dendrite shape selection.

In what follows, we shall need one of the features of equation (8.6), namely

Condition (8.9) appears if we substitute the simplest solution, u(r)=const, of equation (8.2) into which includes the whole contour ζ_i of the domain Ω (figure 4), and generalizes the boundary integral (8.4). Such substitution is correct because the function h(r−r′) depends not on the form of the solution but on the shape of the interface ζ_i.

The kernel (8.6) has the constant α(r′), which is obtained by integration of the delta function along the phase interface,

where r′=ζ_i(r) is the equation of the complete interface and Ω is the liquid region enclosed by the boundary ζ_i, i.e. the region enclosed by the solid/liquid interface ζ_SL and three infinite lines defined by , and (figure 4). For instance, if the interface is a circle with radius R₀, then r′=R₀ is the polar coordinate (scalar) that gives The constant α(r′) appears due to the fact that the argument of the delta function has a null value not in one point (for which the constant α(r′)=1) but on some manifold. In our case, this manifold is represented by the one-dimensional curve (given by the interface) as well as at infinite points.

In obtaining equation (8.4), we have used the transition of the integral with infinite limits for the entire domain ζ_i to the integral on the interface ζ_SL. Indeed, let us assume that the interface ζ_SL propagates in the direction of the z-axis (figure 4). Then, the complete interface ζ_i consists of the phase interface ζ_SL, which represents the solid/liquid interface of a crystal, and three infinite lines defined by , and . In this case, the integral in equation (8.9) is the sum of the integrals over these lines and over the phase interface,

In this sum, the first and third integrals are equal to zero because the Bessel functions are exponentially damping at infinity. Hence, the integral kernel h(r−r′) has the following feature: which will be used in the numerical scheme for obtaining a crystal shape.

(b) Numerical method and modelling algorithm

From the computational viewpoint, the boundary integral method is related to the group of front tracking methods. Even though our equation includes local non-equilibrium terms due to W=V/V_D, it has the form of the equation solved numerically by Saito et al. [31], whose numerical scheme we use to demonstrate the abilities to model dendrites using the boundary integral.

To construct the numerical scheme, let us consider the interface in the form of a polygonal line having 2N+1 points r_j with symmetry relative to the z-axis. Every segment s_j=r_j+1−r_j of the polygonal line is parametrized by λ as

The concentration u and its normal gradient q=∂u/∂n are linearly interpolated between segments at the interface,

Now, the integral equation (8.4) can be approximated by the system of linear algebraic equations, where the matrix elements G_ij and H_ij are computed by the following expressions:

To obtain the diagonal elements H_ii, we have to use equation (8.14), from which one obtains

Using the Gauss method, the integrals of equation (8.18) are approximated by the following finite sum: Weight coefficients w_i and abscissae x_i can be found in the tables in [66]. With computations of the matrix elements having i=j and i=j+1, one should take into account that integrands include the Bessel functions with a singularity in the zero point. For example, For the corresponding matrix elements, we use the Gaussian method of functional integration with the logarithmic singularity The weight coefficients w_i and the abscissae of equation (8.22) are also given in [66].

Using the numerical scheme (8.15)–(8.22), one can model the formation of the steady dendritic pattern using the following algorithm:

• — the solid/liquid interface is chosen initially as a polygonal line that passes through the points r_−N,…,0,…,r_N;

• — the initial value of the normal velocity for these points is W_i;

• — the matrix elements G_ij and H_ij are calculated by equation (8.18);

• — the concentrations in the nodes of the polygonal line, , are calculated using the undercooling balance (9.2), which takes an isothermal approximation into account;

• — in the solution of the system of linear equations (8.17), the values of a normal concentration gradient, q_i, are obtained;

• — the mass balance (8.3) allows us to compute the new values of the velocity W_i of the interface and to choose the new value of the velocity W₀ for the moving reference frame;

• — moving every node in the normal direction on the small distance W_iδt, one obtains the next position of the nodes and entire polygonal line, i.e. the interface;

• — at every point r_i, the curvature is calculated as being inversely proportional to the radius of a circle which passes through this and two neighbouring points; the normal to the interface is chosen along the radius of this circle.

A convergence to the expected approximated solution of equation (8.4) is a special problem for this investigation. Here, it is merely expected that the above-described algorithm and computational scheme provide the convergence to the quasi-stationary shape of a pattern with a constant interface velocity. The described algorithm does not model a process of growth; rather, it just represents iterations to the quasi-stationary dendritic configuration. In this iterative process, the convergence essentially depends on the initial configuration of the interface, the parameter δt and the distance δr between the nodes of the polygonal line.¹ Indeed, δt and δr depend on the length h_D=D_C/V_D, the diffusion length ℓ and the capillary length d₀,

The diffusion length ℓ defines the thickness of the concentration profile, the concentration gradient and the characteristic scale of the steady pattern. Therefore, to adequately predict the dendritic pattern, one has to use the inequality δr≪ℓ. Another limitation comes from the Gibbs–Thomson condition according to which the capillary length defines an action of surface energy for δr>d₀. Thus, the distance between nodes is limited by the double inequality d₀<δr≪ℓ.

9. Pattern selection

The numerical scheme and algorithm formulated in §8 for the boundary integral calculation allow us to study the influence of material parameters and the driving force (undercooling, supersaturation) on the selected quasi-stationary pattern that should appear after the pattern’s temporal evolution. In this work, we demonstrate how the shape of the dendrite depends on the anisotropy of the surface energy and the kinetics of atomic attachment. In addition to this, we shall demonstrate an instability of the Ivantsov analytical solution [40,46] and show that undercooling may be high enough to provide a transition from the non-planar cellular–dendritic front.

The two-dimensional computational domain, in which we simulate convergence to the quasi-stationary dendrite growth (§8), is measured and shown below in units of h_D=D_C/V_D. For the modelling of dendrite formation under constant temperature, we have chosen a Cu-Ni binary alloy with the material parameters given in table 2. In search of the steady-state shape, we quantitatively evaluate the dendrite tip radius R(z)=ρ/2 and dendrite tip velocity W(z)=V (z)/V_D depending on the dimensionless coordinate z measured in units h_D=D_C/V_D.

(a) Selection of dendritic shape

Figure 5 shows the influence of anisotropy on the shape of the dendrite. It can be clearly seen that with the decrease in anisotropy of the surface energy the dendrite becomes thicker and the secondary branching comes closer to the dendrite tip (figure 5a–c). And conversely, the amplitude of the side-branches decreases with the increase in the capillary anisotropy α_d. This effect is known from the experiments of Glicksman and co-workers [93,94].

• Download figure
• Open in new tab
• Download PowerPoint

With the lowest anisotropy (figure 5d), the tip radius oscillates in space, exhibiting a non-constant value. The stationary shape, in this case, does not exist even if the kinetic anisotropy is small but not zero. This example demonstrates that, with the small anisotropy of the surface energy and the relatively small undercooling, the kinetic anisotropy does not compensate the beginning of tip splitting.

As a whole, the steady-state dendrites shown in figure 5 differ significantly from the smooth parabolic solutions of Ivantsov. They exhibit branched patterns, which should appear due to noise imposed along the directions with a maximal value of the anisotropy [95,96]—in our case, this is the initially established parabolic contour. Because we did not use subroutines which smooth the computational noise, one can suppose that, in the present calculations, noise may come from the grid’s discretization and summarizing irregularities arising from the competiting side-branches [31]. With the highly discretizating grid, one can expect that the solution will be smoothed when needle-like growing forms are obtained. In this case, reproducing side-branching will be possible by introducing the fluctuation via thermodynamic noise which can be taken into account by the additional term to the transport flux (thermal or chemical diffusion flux).

(b) Instability of the Ivantsov solution, tip splitting and fractal formation

Figure 6a clearly demonstrates that anisotropy of the surface energy provides a constant velocity and steady tip radius. This result completely confirms one of the main outcomes of the solvability and boundary integral theories developed for the selection criterion and shows that the isotropic Ivanstov solution for a needle-like dendrite [70,71] is unsteady [40,46]. However, the value of the stiffness of the surface energy α_d should be large enough to preserve the steadily growing dendrite tip. Figure 5d shows that, with a very small value of stiffness, the tip of the dendrite has not found its steady form (within the same size of computational domain as in figure 5a–c) and side-branching has already occurred in the vicinity of the tip. At exactly zero value of the anisotropy, the dendrite tip cannot find and select its constant velocity and the tip radius (figure 6b). We observe that the tip splits with branching that, quantitatively, results in spatially non-monotonic behaviour of the velocity W and tip radius R (figure 7). This non-monotonic behaviour of the velocity with the observed tip splitting exists in the formation of six- or threefold symmetrical patterns induced in ferroelectrics [23] as well as in the formation of fractals [97,98].

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(c) Morphological transition at high growth rate

Another effect, which may be demonstrated in applications of the boundary integral theory, is well known as a transition from planar to non-planar growth morphology and vice versa [99]. Here, we demonstrate a transition from the cellular–dendritic pattern to the planar front at high growth rate, which was observed experimentally [83,84,100], predicted theoretically [101] and shown in computational simulations [85–87,102]. We qualitatively and quantitatively describe the transition occurring in isothermal rapid solidification of the binary Cu-Ni alloy melt; see the material parameters in table 2 (which are similar to those in [103], where the transition was modelled using mushy zone equations with a cellular automata algorithm).

The results of the present numerical simulations predict cellular–dendritic patterns (figure 8a) which disappear and transform to the planar front (figure 8b) with an increase in undercooling ΔT=T_QΔ. As in natural and numerical experiments [83–87,100,102,103], cellular–dendritic structures gradually disappear with the appearance of planar-grained forms or macroscopically smooth planar interfaces, forming a ‘structureless’ example of a pattern.

• Download figure
• Open in new tab
• Download PowerPoint

The transition to a planar front in a rapidly solidifying isothermal sample occurs if the phase interface reaches the velocity of absolute chemical stability [104,105]. The two first curves for V and R shown in figure 9 are given by the model predictions from [106] applied to the dendritic solidification of the Cu-Ni melt under constant temperature. With the increase in undercooling, the tip velocity V and tip radius R gradually increase and decrease, respectively. The model [106] predicts that, in close proximity to critical undercooling ΔT_CR, the radius R has a minimum and goes to infinity, which shows the transition to the planar interface (figure 9b). This occurs due to the drastic decrease in the solute boundary layer up to the value of its thickness comparable to the capillary length [104,105]—the scale which stabilizes the planar front. At ΔT=ΔT_CR, the velocity degenerates into a straight line and merges with the velocity

for absolute chemical stability (curve 3 in figure 9a). The velocity W for ΔT≥ΔT_CR is calculated as [107] where the non-equilibrium solute redistribution function k_v and liquidus slope m_v are given by equations (3.3) and (3.5), respectively. The total undercooling ΔT at the planar interface is given by with the solutal undercooling, the undercooling due to the change in the equilibrium liquidus slope m to the non-equilibrium liquidus slope m_v, and ΔT_K=V/μ_k the undercooling which is necessary to attach atoms to the interface. By ΔT_CR, we denote the critical undercooling above which the steady-state growth of the planar interface becomes possible.

• Download figure
• Open in new tab
• Download PowerPoint

Calculations by equations (9.1) and (9.2) show that the critical undercooling for the beginning of the planar interface growth is ΔT_CR≈38 K. From the results of the numerical simulation (figure 8), one can see that the transition from the non-planar to the planar interface proceeds; namely, in the range of ΔT=35−40 K, which includes this predicted value of the undercooling ΔT_CR. Therefore, the present modelling based on the boundary integral theory gives the quantitatively correct interface velocity, as predicted by the sharp interface model [106] based on the solvability criterion for the selected growth mode of the dendrite.

10. Conclusion

The review paper presents analytical and numerical solutions of a Stefan-type thermo-diffusion problem formulated in terms of a single integro-differential equation for the interface function. The previously known results for parabolic heat and mass transfer are detailed and generalized for rapid phase transition processes accounting for the hyperbolic mass transfer for solute diffusion. The generalized integro-differential equation for the transient interface function contains exhaustive information concerning different peculiarities of interface propagation as well as the structural transformations and growing shapes met in the phase transition processes.

The formulated boundary integral equation contains the previously known growing shapes of dendrites obtained by Ivantsov [70,71] and Horvay & Cahn [79] for local equilibrium steady-state solidification regimes (parabolic equation for heat transfer). In the case of rapid solidification processes, when the transport of particles is described by the hyperbolic equation, the boundary integral equation leads to the steady-state solutions with the inifinite rotational symmetry obtained in [73]. If we are dealing with the finite rotational symmetry of an elliptical paraboloid, the solution for the interface function gives the steady-state solution described in [68] for the case of hyperbolic mass transfer, which includes previously known solutions [70–73,79–81] as limiting cases.

In the limit of large Péclet numbers, and , the boundary integral gives semi-spherical solutions for the solid/liquid interface in the vicinity of the dendrite tip. These solutions predict a globular form deformed in the direction of anisotropy strength or a smooth globe in the isotropic case. A qualitative similarity between the boundary integral solution describing globular tips of dendrites and the results of phase-field modelling of globular crystals has been shown for anisotropic and isotropic shapes growing at high growth rates.

In the numerical modelling, we have illustrated the flexibility of the boundary integral theory in reproducing various effects in dendritic growth:

• — selection of the dendrite tip velocity and tip radius;

• — morphological instability of a smooth parabolic ‘crystal/liquid’ interface with secondary branching propagation along the dendrite surface;

• — influence of the interfacial anisotropy on dendrite tip stability with regard to its splitting and the appearance of fractal shapes of crystals;

• — transition from cellular–dendritic structures to the planar front at high growth velocity.

Data accessibility

This article has no additional data.

Authors' contributions

All authors contributed equally to the present review paper.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the Russian Science Foundation (grant no. 16-11-10095) and the German Space Center Space Management under contract no. 50WM1541.

Footnotes

One contribution of 16 to a theme issue ‘From atomistic interfaces to dendritic patterns’.