MODELLING OF FREE BOUNDARY PROBLEMS FOR PHASE CHANGE WITH DIFFUSE INTERFACES

We present a continuum thermodynamical framework for simulating multiphase Stefan problem. For alloy solidification, which is marked by a diffuse interface called the mushy zone, we present a phase filed like formalism which comprises a set of macroscopic conservation equations with an order parameter which can account for the solid, liquid, and the mushy zones with the help of a phase function defined on the basis of the liquid fraction, the Gibbs relation, and the phase diagram with local approximations. Using the above formalism for alloy solidification, the width of the diffuse interface (mushy zone) was computed rather accurately for iron-carbon and ammonium chloride-water binary alloys and validated against experimental data from literature.


Introduction
Free boundary problems (FBPs) involving diffusional and convective transport of heat and mass occur in many engineering problems of practical interest.Some examples of this class of problems are evolution of crystal structure from melt, pulsed laser irradiation, wave interactions on ocean surface, Hele-Shaw cells for pattern formation, seepage of liquid through porous media, and so forth.Cryer [7] has defined the free boundary problems as elliptic and steady-state problems involving dynamically evolving interface as opposed to moving boundary problems which are parabolic and time-dependent.A more exact definition has been put forth by Fasano and Primicerio (see [8,9]).The free boundary problems (FBPs) are characterised by domains which have no known a priori laws for evolution in time.This lack of knowledge regarding the evolution of the phase interface necessitates imposition of an additional condition on the function representing the evolving interface, depending on the physics of the problem.In case of solid-liquid phase change, this additional condition is the energy balance at the phase interface.The mathematical statement of FBP for phase change under thermal diffusion in one dimension is given as follows: find u = u(x,t) and x = Γ(t) for 0 < t < τ such that it satisfies the following conditions: u − u xx = 0, 0 < x < Γ(t), 0 < t < τ, Γ(0) = b, u(x,0) = h(x), 0 ≤ x ≤ b, u Γ(t),t = 0, 0 < t < τ, u x Γ(t),t = −Γ(t), 0 < t < τ, u(0,t) = f (t) or u x (0,t) = f (t), 0 < t < τ. (1.1) The penultimate condition is the energy balance condition at the free boundary and is known as the Stefan condition named after J. Stefan who pioneered the study of FBP pertaining to melting of polar ice caps or conversely freezing of seawater in polar ice caps.
The original formulation of Stefan treated the polar ocean as a column of pure water at bulk equilibrium melting temperature T M adjacent to an ice layer at bulk temperature T b < T M and defined the Stefan condition of latent heat evolution and conduction away from the ice-water interface in a one-dimensional sense [5].FBPs can be further classified as explicit or implicit depending on the explicit or implicit specification of the speed of temporal evolution of the phase interface.The classical Stefan problem falls in the former category.
However, the problems of real engineering and technological interest involving crystal growth or metal casting are not only multidimensional but also multicomponent, multiphase problems.Various fields and fluxes are known to be coupled at various scales of observation in these nonclassical Stefan problems, which govern the complex evolution of the phase interface.These are, for example, transport of mass, momentum, energy, and species at the macroscopic scale (10 −2 − 1m), the solute partitioning at the phase interface at the mesoscopic scale (10 −4 m), undercooling and branching at the dendrite tip at the microscopic scale (10 −5 m), and atomic attachment kinetics at the nanoscopic scale (10 −9 m).Broadly Stefan problems can be classified into (a) phase change of pure materials with a sharp interface, (b) phase change of alloys with a diffuse interface.The solidification and/or melting of ice-water systems under slow cooling rate and phase change materials (PCMs) are examples of the former class whereas solidification of various ferrous and nonferrous alloys, such as steel, cast iron, aluminium alloys, magnesium alloys, and so forth are the examples of the latter class.
In case of the former class of problems, the transient heat diffusion with a source term based on evolution of latent heat has historically constituted the mathematical statement of the problem.This is based on the assumption that heat diffusion alone is responsible for the phase change phenomena.This assumption holds for phase change in pure materials.For phase change in alloy systems, the basic driving force, namely, undercooling, can be achieved both by thermal and solutal diffusion.Their effects could be cooperative or noncooperative at the phase change interface.Further, the advection of heat and solute to and from the interface helps develop microstructural and morphological evolutions at the local and global levels, for example, microsegregation and macrosegregation of solutes, formation of equiaxed or columnar dendritic crystals, and so forth.This necessitates an appropriate framework for analysis where these various fields and fluxes and their mutual interactions are accounted for.
With the help of a modified mixture continuum model based on a previous model [22], we attempt to address the physics for the phase change of alloys.In the mushy zone, the conservation equations are written for the continuum by taking fluxes and forces in a mixture average sense.Exact or analytical solutions for these nonclassical Stefan problems are limited to very simple geometries (e.g., semi-infinite slabs) and uncoupled physics (e.g., phase change driven only by thermal diffusion) [5,2].In this paper, we briefly present the continuum thermodynamic formalism for simulating multiphase Stefan problems with diffuse interface.This is followed by case studies for alloy solidification where the evolution of the phase front is tracked by numerical simulation and compared with reported experimental data.

Stefan problem with a diffuse interface
For an alloy system undergoing phase change, the mushy zone consists of a network of solids of columnar, equiaxed, or mixed dendritic morphology.The conventional Stefan problem fails to model such a system and the free boundary can become unstable rendering the problem as ill-posed.This needs reformulating the problem with the introduction of a phase function describing the mushy zone.In Atthey [1], a model problem of heat conduction in a slab was considered where melting was produced by distributed heat sources.The mushy zone was characterised by the heat content.It was thought that the mush was formed due to addition of heat which was insufficient for complete melting.The weak solution of this problem was computed which showed the appearance of a finite region where the temperature was equal to the melting temperature.The smoothness of the weak solution was proved by Primicerio [19] by solving a differential system in a classical sense.Lacey and Tayler [16] developed a very simple model for phase change with fine internal structure of some mush which upon averaging yielded the macroscopic formulation of Atthey [1].Ughi [23] considered a similar formulation as that of Atthey [1] and investigated some relevant properties of the solution for a one-dimensional problem in a slab.Visintin [24] proposed a generalisation of the Stefan problem for phase transitions in one-dimensional systems taking account of the nonequilibrium supercooling and superheating effects for sharp and diffuse phase fronts.They proved the existence of at least one solution by means of Faedo-Galerkin approximation procedure and gave some complementary results.Their general formulation yielded the standard one-phase Stefan problem as a limit case.However, the framework of rational thermodynamics and multiphase dynamics has emerged as the basis of some more general models of phase change in general and mushy zone in particular [15,25].
Under the assumptions of relatively slow cooling of the melt, local thermodynamic equilibrium is expected to hold at the interface.The accurate description of the thermodynamics of the phase interface is governed by the phase diagram and the Gibbs relation in conditions close to equilibrium.For binary alloy solidification comprising the solute γ and solvent α, the phase diagram can be represented in terms of the composition of the solute C γ = ζ ± (T) where superscript + represents the solidus, superscript − represents the liquidus curves, and ζ is a continuous function.
Considering an elemental volume dΩ which is small and is characterised by uniform properties, its properties are determined by the magnitude of the function ζ.Without loss of generality, it can be shown that the volume dΩ belongs to a pure solid, pure liquid, or a mushy zone of the continuum as where Ω + , Ω − , and Ω M denote the solid, liquid, and the mushy zones, respectively.The mushy zone is a part of the continuum which is not only a superposition of Ω + and Ω − , but also comprises the chemical species α and γ constituting the binary alloy.The membership of dΩ into various domains of the Ω can be conveniently described in terms of the liquid fraction f L as follows: where is the equilibrium solute partition ratio at temperature T.
The Lever rule assumes infinitely fast diffusion of the solute species from the interface to the bulk in both the liquid and solid phases whereas the Gullivers-Scheil equation assumes zero diffusion in solid phase.Any property ψ (such as the enthalpy, internal energy, chemical potential, etc.) defined for the entire continuum can be written as where • is the jump condition and is defined as follows: Before we write the conservation equations for the mass, momentum, energy, and species transport, it is imperative that we describe the thermodynamics of the phase change process using the equilibrium phase diagram, as given by (2.1) and the Gibbs relation for the liquid and solid phases, given as follows: where s is the specific entropy, i 's are the generalised forces, and χ i 's are the generalised coordinates.For the binary alloy solidification 1 dχ 1 = pdV and 2 dχ 2 = −(g γ − g α )dC γ , where p is the pressure, V is the specific volume, and g ν is the specific Gibbs free energy for species ν.
From the classical thermodynamic definition, the specific internal energy can be defined as a function U = U(T,V ,C γ ) such that the following equation of state results: From the thermodynamic relations, ) ) ) ) where g ν and h ν are the specific Gibbs free energy and the specific enthalpy for the species ν and ω is the difference in specific enthalpies of species γ and α, respectively.Substituting (2.7) to (2.11) in (2.6), the differential specific internal energy can be written, with the help of the Gibbs-Duhem equation and after several algebraic manipulations (refer [15,25] for details), not shown here for paucity of space, as follows: and M), (2.12) where superscript k implies the phase.∆H f is the latent heat of fusion.In pure liquid and solid zones df ± L = 0. Due to the locally isobaric condition, the above equation can be rewritten as where ϑ k and ξ k are the thermal and solutal expansion coefficients of the kth phase.The conservation equations for the continuum representing the solidifying binary alloy is given by the following equations.
(i) Continuity: (2.15) V , V + , and V r are the mixture velocity, solid phase velocity, and the relative velocity of the mixture with respect to the solid phase, respectively.µ eff is the effective viscosity of the fluid phase, K y is the y-component of the permeability tensor, and N is a constant in the Forscheimer term.Similar equations can be written for conservation of the momentum for the x-component and z-component, respectively.
(iii) Energy: (2.16) (iv) Species: (2.17) In (2.16), D is the diffusivity of the alloying species and Θ is the radiation heat source.In (2.16) and (2.17), terms involving δ and β are due to Soret's effect and the Dufour effect, respectively.The term ∇.[(ρD/C α ){∂(g γ − g α )/∂C γ } −1 ] is defined as the mobility and when multiplied by the acceleration due to gravity (g), represents forced diffusion due to gravitational forces.For detailed derivation of the above equations [25] may be referred to.However, some of the terms in the above equations need further explanation.Equation (2.15) is the Navier-Stokes equation wherein the fourth and fifth terms on the right-hand side represent the Darcy and Forscheimer phase interaction terms, respectively.The last term of (2.15) represents the Obverbeck-Boussinesq approximation accounting for thermal and solutal buoyancy components.In absence of the Soret and Dufour effect and for negligible mobility, (2.16) and (2.17) can be conveniently simplified to the more familiar mixture continuum forms derived in [11].Equations (2.14)-(2.17)represent a set of parabolic partial differential equations with strong coupling between thermal, velocity, and concentration fields.In addition, these need to be coupled with an appropriate model of turbulence in case of turbulent flow.These are not amenable to analytical solutions and need numerical solution.Because of the necessity to track the interface, the Stefan problem needs special algorithms as briefly outlined below.

Interface tracking scheme
Stefan problems with diffuse interfaces have been simulated either with (a) multidomain, deforming grid algorithm or (b) single-domain, fixed grid algorithm.In the former approach, separate sets of conservation equations are written for the solid and liquid phases and a heat balance equation is written at the interface which is known as the multiphase Stefan condition.In the single-domain approach, a single set of macroscopic conservation equations holds for all the phases with a phase field or a phase function accounting for the membership of the control volume to a particular phase.The single-domain approach is more suited to the solution of alloy solidification problems.For numerical solution of single-domain models of phase change, the fixed grid strategies are usually preferred as they retain the same Cartesian mesh throughout computation to describe the motion of the interface over time and space.In this category, three major techniques stand out for their popularity and ease of use, namely, the finite difference/finite volume method, the finite element method, and the boundary element method, though, not necessarily in that order.Typically for a finite difference or finite volume discretisation, a mesh point or a control volume surrounding a mesh point (or node) is assigned a value of the phase function (linked to the volume fraction of solid) based on the local thermal field and the thermodynamics of phase change.The propagation of the interface is then captured by the phase function as it evolves over the domain in space and time.However, based on the treatment of the evolution of the latent heat at the interface, several broad classes of fixed grid strategies have been developed in the literature.Poirier and Salcudean [18] have identified eight such categories.In the postiterative method [18] for a pure material, an enthalpy budget is maintained for the finite difference or finite element nodes at which phase change occurs.At these nodes the temperature is periodically set to phase change temperature and an equivalent amount of latent heat is added to the node's enthalpy till it equals the latent heat for the volume associated with that node.The temperature is allowed to fall according to heat diffusion.An extension of the algorithm for mushy zone has been proposed by Salcudean and Abdullah [20] where the nodes falling within the mushy range are set to a temperature where subscript "liq" refers to liqidus and "sol" refers to solidus.∆H is the heat content for the enthalpy budget and ∆H f is the latent heat of fusion.
In apparent heat capacity technique, the release of latent heat at the interface is accounted for by a suitable increase in the heat capacity in the range of temperature pertaining to phase change.For a linear release of latent heat across this temperature range, the apparent heat capacity is given as C pl for T > T liq .Since enthalpy variation is smooth across a phase change interface in comparison with the apparent heat capacity (see Figure 3.1), the enthalpy method [6,10] was developed to exploit this feature such that the governing equation is recast in terms of enthalpy as follows: where the enthalpy H, given as follows, is a continuous function of temperature T except at the melting point of a pure material or at the eutectic temperature for an alloy where a jump in enthalpy occurs (Figure 3.1): The major advantage of the fixed grid approach is that while it tracks the large distortions in the evolving phase front, the grid topology remains fixed and simple.This makes the computation less expensive and requires much less memory resources.In the present paper, a fixed grid control volume finite difference scheme has been employed.A higher-order upwinding formulation has been proposed and embedded within the overall scheme based on the SIMPLER of Patankar [17].The details of the scheme is available in [21] and is excluded from the present paper for the sake of brevity.

Numerical simulation
The above mathematical formalism of Stefan problem based on advection-diffusion of the alloying species coupled with energy and momentum transport along with two equation turbulence model has been employed for numerical simulation of the solidification of two binary alloys (a) hypereutectic steel of nominal composition 0.8 weight percent carbon and (b) hypereutectic ammonium chloride-water binary alloy of nominal composition of 30 weight percent ammonium chloride.The former is a short freezing range alloy of iron and carbon.For case (a), solidification of steel is considered in a slab caster of rectangular cross-section 1.6m × 0.22m.The total slab length simulated is 3m from the meniscus.The molten steel is poured through a nozzle 0.03m in diameter from the top of the caster.The length of the mould is 0.5m and the caster speed is 0.01m/s.The operating and thermophysical conditions of case (a) are given in Table 4.1.The latter system, that is, NH 4 Cl-H 2 O system, is used for experimental studies popularly for its optical transparency which makes it amenable to visual observation of the evolution of the solidification interface and the hydrodynamics.The solidification was done in a three-dimensional cavity of rectangular cross-section of L = 36mm and H = 144mm (aspect ratio, A = 4.0) with a thickness of 2W = 200mm.The computational domain was a rectangular section in two dimensions of 36mm × 144mm.The simulation runs were conducted for various operating conditions for eutectic, hypereutectic, and hypoeutectic alloy compositions.These are reported in detail in [21].We present here a case study for a simulation run where eutectic and hypereutectic alloy of nominal composition of 20 and 30 weight percent ammonium chloride, respectively, is allowed to solidify from an initial temperature of 40 • C when the left wall of the cavity was maintained at −30 • C and the right wall was insulated.The thermophysical data are reported in Table 4.2.
With respect to a computational domain discretised with an optimum nonuniform mesh size of 40 × 22 × 100 for the steel slab and 46 × 74 for the ammonium chloridewater binary alloy, the simulation runs were conducted for steady-state operation in case of casting of steel slab and for various times from the start of the solidification for ammonium chloride-water binary, respectively.The growth of solid steel shell and the mushy zone under steady-state condition is given in Figure 4.1a.As steel is a short freezing range alloy of iron and carbon, the mushy zone width is very narrow.The computed shell thickness is in general agreement with industrially observed slabs cast with the above composition [13,14].Figure 4.1b depicts the computed contours of the solid shell (white zone) and a very negligible thin band of mushy zone in shades of gray near the interface of liquid NH 4 Cl-H 2 O melt (dark zone) and the solid shell.This is expected for an alloy of near eutectic composition which is characterised by an absence of the mushy zone.This goes to show that the diffuse interface Eulerian fixed grid algorithm is adequate in handling the evolution of nearly sharp interface with a simple computational algorithm in comparison with the more complex Lagrangian interface tracking algorithms.For the hypereutectic alloy composition, a wide mushy zone is expected to form because of the large freezing range indicated from the alloy phase diagram.This is observed in the photographic images reported for NH 4 Cl-H 2 O binary in [3]. Figure 4.1c depicts the image of the solidus front (bright shade near the cold wall on the left) and liquidus(mushy zone) front (dark shades with rough contours) for the same alloy composition after 22 minutes from the start of solidification.The computed contours of the solid shell(solidus) and mushy zone(liquidus) are depicted in Figure 4.1d.A black line is drawn as an isoconcentration contour at solid fraction of 0.1 (approximately) marking the boundary between the mushy zone front and the liquid alloy.The match between the experimentally observed and the computed solid fraction is reasonably good with absolute average error of 8%; much better than that computed by Christensen and Incropera themselves [4] which underpredicted the mushy zone width by 10-36% with absolute average error of 26%. Figure 4.2 depicts the comparison between the experimentally estimated mushy zone front and that computed in this paper as well as by Christensen and Incropera [4] for the same alloy composition with the right wall at 40 • C.
To highlight the efficacy of the numerical technique, validation studies have been done with simpler test cases.In this paper, we present one such test case which pertains to solidification of ice in a transparent cubical cavity made of plexiglass with each side 0.038m long in an experimental study reported by Giangi et al. [11,12].The hot wall (on the left) was maintained at 10 • C and the cold wall (on the right) was initially maintained at 0 • C. The flow was allowed to attain steady state till 120seconds after the onset of convection.The thermal and velocity fields in the water were measured using thermochromic liquid crystal (TLC) tracers, digital particle velocimetry, and thermometry.The most interesting  The study shows that the flow is relatively insensitive to grid refinement near the hot wall (X = 0.066), but has a strong sensitivity near the cold wall (X = 0.9) due to the density inversion.A uniform grid of 42 × 42 was found acceptable as a compromise between accuracy and computational cost.

Conclusions and future directions
In this paper, a thermodynamically consistent continuum mixture formalism has been presented as a suitable framework for modelling alloy solidification phenomena.Numerical simulations of solidification phenomena in binary alloys of iron-carbon (steel) and ammonium chloride-water have been performed with a fixed grid algorithm and have been validated, for the latter, with experimental data of Christensen and Incropera [3].From the close match between the mushy zone width computed in this paper and that obtained experimentally [3], the present algorithm is proven as effective and an improvement over the previous model of Christensen and Incropera [4].The efficacy of the algorithm has been demonstrated with a simpler case study on convective flow during freezing of ice [12].Further simulation runs are being conducted with different sets of alloy solidification data for more rigorous validation of the same and will be reported in the future.

Figure 4 .
3 shows an excellent match between (a) the velocity streamlines experimentally observed and (b) numerically computed in this study.Grid refinement studies have been carried out to test the efficacy of the numerical schemes.This is reported in Figure4.4.The velocities plotted in Figure4.4 have been nondimensionalised using a characteristic velocity of α/H = 3.545 × 10 −6 m/s.

Figure 4 . 1 .
Figure 4.1.Solid shell and mushy zone width in steel and NH 4 Cl-H 2 O binary alloy.(a) Computed profile for steel.(b) Computed profile for eutectic NH 4 Cl-H 2 O. (c) Experimentally observed profile in 30 weight percent NH 4 Cl-H 2 O [3].(d) Computed profile for 30 weight percent NH 4 Cl-H 2 O. (All computed profiles pertain to the present work.)

Figure 4 . 2 .
Figure 4.2.Mushy zone width in NH 4 Cl-H 2 O binary for 30 weight percent nominal composition and right wall at 40 • C.

Figure 4 . 3 .Figure 4 . 4 .
Figure 4.3.(a) Black and white intensity images from TLC tracers reported in [12].(b) Calculated streamlines for the flow obtained in the present paper.

Table 4 .
1. Operating and thermophysical parameters for steel slab casting.
and challenging aspect of simulating this problem is to accurately model the structure and location of the two opposing buoyant convection currents due to density inversion of water at 3.98 • C.