Mathematical models and numerical solutions of liquid-solid and solid-liquid phase change

1.1. Introduction

The phase change phenomena in which the matter transitions and transforms from one state to another is of significant academic and industrial importance. Solid-liquid or liquid-solid phase transitions and their numerical simulation have been a subject of research and investigation for over a century. There are many sources of difficulties in the numerical simulation of phase change phenomena. Phase transition physics and its mathematical modeling is quite complex due to the fact that this phenomenon creates a transition region, a mixture of solid and liquid phases, in which the phase change occurs resulting in complex changes in transport properties such as density, specific heat, conductivity and the latent heat of fusion that are dependent on temperature. During evolution the phase transition region propagates in spatial directions, i.e. its location changes as the time elapses. Idealized physics of phase change, in which jumps in the transport properties are often assumed, results in singular interfaces. As a consequence the mathematical models describing such evolutions result in initial value problems that contain singularities at the interfaces. When solving such non-linear initial value problems, one must assume existence of the interface. Numerical simulation of the propagation of such fronts during evolution also presents many difficulties that cannot be resolved satisfactorily. Major shortcomings of this approach are that formation of the phase transition front cannot be simulated. Secondly, singular nature of the front is obviously not possible to simulate numerically. In the second approach of phase transition physics and its mathematical modeling, one assumes that the phase transition region is of finite width, i.e. the phase transition occurs over a finite but small temperature range in which the transport properties such as density, specific heat, conductivity and latent heat are function of temperature and vary in a continuous and differentiable matter between the two states. Thus, the phase transition region is of finite width in temperature that propagates as time elapses. This approach is more realistic and more appealing from the point of view of numerical simulations of the resulting IVPs from the mathematical models as it avoids singularities present in the first approach. The phasefield approach utilizes this concept. A major source of difficulty in this approach is the physics of the transition region, often referred to as ‘mushy region’, that consists of liquid-solid mixture in varying volume fractions as one advances from one state to the other. Adequate mathematical modeling of the physics in the transition region may require use of mixture theory [1–3] or some similar approach, based on thermodynamic principles of continuum mechanics. Conservation of mass, balance of momenta, first law of thermodynamics and the constitutive theories for stress tensor and heat vector based on the second law of thermodynamics must all be reformulated assuming thermodynamic equilibrium in the transition region. This approach of mathematical modeling of the transition region has not been explored in the published literature (to

of the methods of approximation that are utilized to obtain numerical solutions of the initial value problems describing evolution. It is now well established in computational mathematics that methods of approximation such as finite difference, finite volume and finite element methods based on Galerkin Method (GM), PetrovGalerkin method (PGM), weighted residual method (WRM), and Galerkin method with weak form (GM/WF) used in context with space-time decoupled or space-time coupled methodologies are inadequate for simulating time accurate evolutions of the non-linear IVPs describing phase change processes [4–9]. Thus, in order to address numerical solutions of phase transition processes, in our view a simple strategy would be to: (i) Decide on a mathematical model with desired, limited physics. (ii) Employ a method of approximation that does not disturb the physics in the computational process, results in unconditionally stable computations and has inherent (built in) mechanism of the measure of error in the computed solution without the knowledge of theoretical solution as such solutions may not be obtainable for the problem of interest. The work presented in this thesis follows this approach. In the following we present literature review on mathematical modeling and methods of approximation for obtaining numerical solutions of the IVPs resulting from the mathematical models. This is followed by the scope of work undertaken in this paper.

1.2. Literature Review

In this section we present some literature related to liquid-solid and solid-liquid phase transition phenomena. We group the literature review in two major categories: mathematical models and methods of approximation for obtaining numerical solutions of the initial value problems resulting from the mathematical models. 1.2.1

A large majority of published work on the mathematical models for phase change processes consider Lagrangian description only, with further assumptions of zero velocity field, i.e. no flow and free boundaries i.e. the medium undergoing phase change to be stress free. We first present literature review and a discussion of commonly used mathematical modeling methodologies in Lagrangian description based on the assumptions stated above. With the assumptions of no flow and stress free medium, the mathematical model of the phase change process is invariant of the type of description and reduces to the energy equation. In the published works there are three commonly used approaches: sharp-interface models, enthalpy models and phase field models. In the mathematical models derived using sharp-interface the liquid and solid phases are assumed to be separated by a hypothetically and infinitely thin curve or surface called sharp interface or phase. The transport properties such as density, specific heat and conductivity are assumed to experience a jump at the interface. 63

The latent heat of fusion is assumed to be instantaneously released or absorbed at the interface. This of course results in step (sharp) change in the transport properties and latent heat of fusion at the interface, hence the name sharp-interface models. The mathematical models for liquid and solid phases are derived individually. At the interface, the energy balance provides an additional relation (equation) that is used to determine the movement of the interface. The sharp-interface models are also called Stefan models, first derived by J. Stefan [10] to study freezing of ground. The derivation of this model is presented in Section 2. The proof of existence and uniqueness of the classical solution of the Stefan mathematical model has been given by Rubinstein [11] in 1947. An analytical solution for temperature for one dimensional Stefan problem has been presented in reference [12]. The sharp-interface models have three major shortcomings: (i) Assumption of sharp-interface leads to mathematical model in which the initial value problem contains singularity at the interface. (ii) When obtaining solutions of the initial value problems based on sharp-interface assumption, the location of the interface is required a priori. That is sharp-interface models are unable to simulate initiation of the interface or front. (iii) Movement of the interface i.e. spatial location during evolution requires use of what are called front tracking methods. Some mathematical models for phase change processes are called enthalpy models. In these models the energy equation is recast in terms of enthalpy and temperature with an additional equation describing enthalpy. Both enthalpy and temperature are retained as dependent variables in the mathematical model. Computations of the numerical solution of the resulting initial value problem are performed on a fixed discretization. This approach eliminates energy balance equation at the interface used in the sharpinterface models. These mathematical models have been derived using different approaches [13–15]. Enthalpy model is also presented in Section 2. These models generally introduce a finite phase transition region (over a small temperature change) called mushy region between the liquid and the solid phases. The transport properties are assumed to vary in some manner from one phase to the other phase. The concept of liquid or solid fraction is generally introduced to account for the fact that the mushy region is a mixture of solid and liquid phases. Due to the assumption of the mushy region separating the solid and the liquid phases, sharpinterface and the problems associated with it are avoided in this approach. Another category of mathematical models are called phase field models. These mathematical models are based on the work of Cahn and Hilliard [4]. In this approach the solid and liquid phases are also assumed to be separated by a finite width (in temperature) transition region in which the transport properties are assumed to vary with temperature between the two states. Landau-Ginzburg [5] theory of phase transition is used to derive the mathematical model. The basic foundation of the method lies in standard mean theories of critical phenomena based on free energy functional. Thus, the method relies on specification of free energy density functional which is the main driving force for the movement of the phase transition region. Details of phase field mathematical model in R1 are

presented in section 2. The method shows good agreement with the Stefan problem in R1 . While the phase field models eliminate the sharp-interfaces and their tracking, the main disadvantages of this approach are: (i) It requires a priori knowledge of the free energy density functional for the application at hand. (ii)The mathematical model is incapable of simulating the initiation or formation of the solid-liquid interface, hence the liquid-solid phases and the transition region must be defined as initial conditions. This limitation is due to specific nature of the free energy function (generally a double well potential, see section 2). However, if a liquid-solid interface is specified as initial condition, then the phase field models are quite effective in simulating the movement of the front during evolution. In most applications of interest, simulation of initiation of the transition region i.e. solid-liquid interface is essential as it may not be possible to know its location and the precise conditions under which it initiates a priori. These limitations have resulted in lack of wide spread use of these mathematical models in practical applications. When the assumptions of stress free media and zero velocity are not valid (as in case of fluid flow), the mathematical models discussed above are not applicable. In such cases Eulerian description is necessary for the fluid while Lagrangian description is essential for the solid region. The mathematical model in this case consists of conservation of mass, balance of momenta, first law of thermodynamics and constitutive theory for stress tensor and heat vector based on the second law of thermodynamics for each of the two phases (i.e. liquid and solid) as well as the transition region. The published works on these mathematical models are rather sketchy, the models are not based on rigorous derivation and in most cases are aimed at solving a specific problem as opposed to developing a general infrastructure that addresses totality of a large group of applications. We present some account of the published works in the following. In almost all cases the fluid is treated as Newtonian fluid. In some cases [16] the fluid is also considered inviscid. Sharp-interface models generally force (set) the relative movement of the material particles to be zero in the solid phase [17, 18]. In case of enthalpy and phase field models the constitutive theory for the transition region is still unclear and published works in many instances are conflicting. There are three main ideas that are commonly found in the majority of the published works on mathematical models derived using Eulerian description. In the first approach both the liquid and the solid phases are assumed to be Newtonian fluids. The viscosity in the solid phase is artificially increased to a very high value and is assumed to vary along the interface between the two states in order to approximate no velocity condition in the solid phase [19]. In the second approach a varying interfacial force is employed such that it satisfies the no velocity condition in the solid phase [20]. The third approach assumes that the solid particles in the transition region form a porous medium through which the fluid flows. Voller and Cross [15] use Darcy model for flow in porous media in which the velocity field is assumed to be proportional to the pressure gradient in order to compare their results with variable viscosity model. Beckermann [21] assumed the average stress to be proportional to the gradient of

superficial liquid viscosity in the porous media. There are other approaches [22] that utilize these three basic ideas in some manner or the other. In most cases, solid phase behavior is neglected by setting the velocity to zero. In general, our conclusion is that published phase change models that account for nonzero stress and velocity fields are crude, ad hoc and are aimed to obtain some numerical solutions for specific applications. A general theory of mathematical modeling based on thermodynamic and continuum mechanics principles is not available for phase transition modeling to our knowledge. 1.2.2

Regardless of the type of mathematical model, the resulting mathematical models for phase change phenomena are non-linear partial differential equations in dependent variables, space coordinates and time, hence they are non-linear initial value problems. If we incorporate realistic physics of phase transition, the mathematical models become complex enough not to permit determination of theoretical solution, hence numerical solutions of these IVPs based on methods of approximation are necessary. The methods of approximation for IVPs can be classified in two broad categories [6–9] : space-time decoupled methods and space-time coupled methods. In space-time decoupled methods, for an instant of time, the spatial discretization is performed by assuming the time derivatives to be constant. This approach reduces the original PDEs in space and time to ODEs in time which are then integrated using explicit or implicit time integration methods to obtain evolution. Almost all finite difference, finite volume and finite element methods (based on GM/WF) used currently [7] for initial value problems fall into this category. The assumption of constant time derivatives necessitates extremely small time increments during the integration of ODEs in time. The issues of stability, accuracy and lack of time accuracy of evolution are all well known in the space-time decoupled approaches. Majority of the currently used methods of approximation for phase change processes fall into this category. The non-concurrent treatment in space and time in space-time decoupled methods is contrary to the physics in which all dependent variables exhibit simultaneous dependence on space coordinates and time. In a large majority of published works on phase change processes, often the distinction between the mathematical models and the computational approaches is not clear either i.e. elements of the methods of approximation are often introduced during the development of the mathematical models. As a consequence, it is difficult to determine if the non-satisfactory numerical solutions are a consequence of the methods of approximation used or the deficiencies in the mathematical models. The space-time coupled methods on the other hand maintain simultaneous dependence of the dependent variables on space coordinates and time [6, 8, 9]. In these methods the discretizations in space and time are concurrent as required by the IVPs. These methods are far superior to the space-time decoupled methods in terms of mathematical rigor as well as accuracy. Whether to choose space-time finite difference, finite volume or finite element method

depends upon the mathematical nature of the space-time differential operator and whether the computational strategy under consideration will yield unconditionally stable computations, will permit error assessment, and will yield time accurate evolution upon convergence.

1.3. Scope of Work

The work presented here considers development of mathematical models and their numerical solutions for solid-liquid and liquidsolid phase transition of homogeneous, isotropic, and incompressible matter. In the phase transition region [Ts , Tl ] the matter is assumed to be homogeneous and isotropic and the transport properties are assumed to be continuous and differentiable with their respective values at the solid and liquid states. Three groups of mathematical models are considered for phase transition initial value problems. Numerical studies are presented using the mathematical models groups one and three. The first group of mathematical models are based on the assumptions of stress free media and zero velocity in all phases. With these assumptions the mathematical models in Lagrangian and Eulerian descriptions are identical. We consider these mathematical models in R1 and R2 . The mathematical models in this case consist of the energy equation and heat flux(es), a system of first order nonlinear PDEs in temperature and heat flux(es). By substituting heat flux(es) into the energy equation the mathematical model can be reduced to a single non-linear diffusion equation in temperature. In the derivation of the energy equation the specific total energy is expressed in terms of storage and latent heat of fusion. The Fourier heat conduction law is assumed to hold. In the solid and liquid phases the transport properties (ρ, c p , k, L f ) are assumed to be constant. In the transition region the solid-liquid mixture is assumed to be isotropic and homogeneous. The transport properties are assumed to vary in a continuous and differentiable manner, described by a third or a fifth degree polynomial with continuous temperature derivatives at the boundaries of the transition region between the solid and liquid phases. With this approach the phase change process is a smooth process in which the transition region provides the smooth interface. We remark that if we assume both phases to be incompressible, then a change in density during phase change must be accompanied by a change in volume. In the present work we consider phase change studies in R1 and R2 assuming (i) the density ρ to be constant during the phase transition and (ii) the density to be a function of temperature i.e. variable with continuous and differentiable distribution between the states. Additionally, the influence of temperature dependent density in the transition region on the speed of propagation of the transition region is also investigated. Mathematical models and numerical studies are presented in R1 and R2 for solid-liquid and liquid-solid phase change when stress field and velocity field are zero. In the second group of mathematical models stress and velocity fields are considered to be nonzero. In this case the mathematical models change drastically compared to the first group of models. This is due to the fact that in solid regions Lagrangian 65

2.1. First group of mathematical models for

phase change based on zero stress and velocity fields and free boundaries These mathematical models constitute the first group of mathematical models. When the media are stress free, the velocity field is zero, and the boundaries are free the mathematical model for phase change reduces to linear or nonlinear diffusion equation regardless of the choice of dependent variables. In the published works there is a lot of confusion in the presentations of these models regarding the choice of conflicting notations, representation of physics, and even consistency of derivations. These models are generally classified as sharp-interface models, enthalpy models, phase field models, smooth-interface models, etc. We show that the energy equation resulting from the first law of thermodynamics is the same in all of these models. What differs is (i) the choice of dependent variable(s) and (ii) the manner in which the phase transition physics is incorporated. We present two basic forms of the energy equation that are used in the mathematical models mentioned above. Overbar on quantities indicates that the description is Eulerian with transport. Energy Equation Following [23] for a compressive and dissipative medium, we can derive the following energy equation from the first law of thermodynamics in Eulerian description with transport when the stress field and the velocity are not zero. Assuming sources and sinks to be absent

The third group of mathematical models are derived based on the assumption that the stress field is constant and the velocity field is zero in the solid region but nonzero in the liquid region. In the transition zone, the stress and the velocities are assumed to make transition from nonzero state in the fluid to constant stress state and zero velocity in the solid phase based on the temperature in the transition zone. These mathematical models permit phase transition studies in the presence of flow, are consistent description based on continuum mechanics, and hence provide correct interaction between the solid and fluid media. Numerical studies are presented in R1 and R2 to demonstrate various features of the mathematical models presented here. Computed solutions in R1 are also compared with sharp-interface theoretical solution.

ρ̄ is density, ē is specific internal energy, q̄q is the heat vector, [σ̄ (0) ] is the contravariant Cauchy stress tensor, and [D̄] is the symmetric part of the velocity gradient tensor, all in the current configuration at time t. Equation (1) can also be written in terms of specific enthalpy h̄. Recall that h̄ = ē +

Computed mathematical solutions reported in the paper are always converged and are independent of mesh size and degree of local approximation. In all cases the integrated sum of squares of the residuals are small (O(10−6 ) or lower), confirming good accuracy of the reported solutions. 66

Consider decomposition of [σ̄ (0) ] into equilibrium stress p̄[I] For this case h = e as obvious from (2) when p = 0. In the and deviatoric stress [d σ̄ (0) ] energy equations (12) and (13) the simplest constitutive theory for heat vector is of course Fourier heat conduction law. σ (0) σ (0) = − p̄II + d σ̄ σ̄ (4) Using (4)  
  
  
tr σ̄ (0) D̄ = −tr p̄ D̄ + tr d σ̄ (0) D̄ ¯ · v̄v + tr  σ̄ (0) D̄
= − p̄∇

in which k is the thermal conductivity for homogeneous isotropic matter. Equations (12) and (14) or (13) and (14) form the basis for phase transition mathematical models in the absence of stress field and velocity field. Various methods published in the literature differ in the manner in which the phase change physics is incorporated in (12) and (13).

(1) First we note that since h = e, the specific enthalpy and the specific energy models are the same. From now onwards, we will use (12) to present further details.

∀(x̄x ,t) ∈ Ωx̄x t = Ωx̄x × Ωt (2) The fundamental issue is the physics for e we wish to consider during the phase change. We consider two possibilities.

  
Dh̄ D p̄ − − p̄∇¯ · v̄v + ∇¯ · q̄q + p̄∇¯ · v̄v − tr d σ̄ (0) D̄ Dt Dt ∀(x̄x ,t) ∈ Ωx̄x t = Ωx̄x × Ωt

∀(x̄x ,t) ∈ Ωx̄x t = Ωx̄x × Ωt (9) D p̄ in (9) is often neglected if compressibility is not significant. Dt

(a) In the first class of mathematical models we assume that the release or absorption of latent heat during phase change occurs at a constant temperature. Referring to figure 1(a) when the temperature in the solid medium reaches Ts with specific internal energy es (point B), the addition of latent heat of fusion L f at constant temperature Ts increases es to el (point C) at which the state of the matter has changed from solid to liquid. In case of freezing we go from the state of the matter at C to B by extracting latent heat of fusion L f at constant temperature Ts .

Equations (1) and (10) are two fundamental forms of the energy equation in specific internal energy ē and specific enthalpy h̄ when the medium is compressible and the stress field and the velocity field are not zero.

In this physics of phase transition the interface between the solid and the liquid phases is sharp (step change), hence the mathematical models for e based on this approach are called ”sharp-interface models.” Step change in e is nonphysical even for the most idealized materials.

Secondly, its numerical simulation poses difficulties due to non-unique behavior of e at temperature Ts . We present details of sharp-interface models in a following section.

When the medium is stress free and the velocity field is zero then ∂ D σ (0) = 0 = D̄ D = and (11) d σ̄ Dt ∂t Furthermore, with these assumptions Eulerian and Lagrangian descriptions are the same, hence the overbar on all quantities can be omitted. Thus, (1) and (10) reduce to ρ

(b) In the second category of mathematical models for e we assume that phase transition from solid to liquid occurs over a finite but small range of temperature [Ts , Tl ] and that e is continuous and differentiable for Ts ≤ T ≤ Tl (figure 1(b)). The range [Ts , Tl ] can be as narrow or as large as desired. The obvious advantage in this approach is that the singular nature of e at Ts (as in figure 1(a)) is completely avoided. This is of immense benefit in numerical computations of evolution of phase change problems. 67

Figure 1: Sharp- and smooth-interface models for specific internal energy

Ωs and Ωxl are solid and liquid spatial domains. Γx (t) = Tx l s Ωx Ωx is the interface between the solid and liquid phases. L f

As described earlier these models for e are based on its behav- is the latent heat of fusion, n is the unit exterior normal from the ior during phase transition shown in figure 1(a). We have some solid phase at the interface, and vn is the scalar normal velocity of alternative forms of the mathematical models. the interface in the direction of n . Subscripts and superscripts s and l stand for solid and liquid phases. When the mathematical model is posed as a system of inteModel (a) gral equations, a complete proof of existence and uniqueness of In this model we consider the classical solution in R1 was given by Rubinstein in 1947 [11]. For the one dimensional case, analytical solutions to some specific ∂e ∇·q = 0 +∇ ∀(xx ,t) ∈ Ωx t = Ωx × Ωt (15) problems are derived in reference [12] for the temperature distribuρ ∂t tion T = T (x,t). When the properties are the same in both phases ∇T q = −k∇ ∀(xx ,t) ∈ Ωx t = Ωx × Ωt (16) (i.e. c ps = c pl = 1, ks = kl = 1, ρs = ρl = 1), one example problem in reference [12] solves for T in the domain x ≥ 0 with initial and e = es + αL f + c p (T − Ts ) (17) boundary conditions:  0 ; e < es   e−e T (0,t) = T0 s ; es ≤ e ≤ es + L f (18) α= T (x, 0) = Θ(x) (23)  Lf  1 ; e > es + L f T (x,t)x→∞ = T∞ Alternatively (16) can be substituted into (15) to obtain Then the solution to the sharp-interface model is given by     ∂e x β ∇ · (k∇ ∇T ) = 0 ρ −∇ ∀(xx ,t) ∈ Ωx t = Ωx × Ωt (19) − erf √ erf ∂t 2 2 t + t0   ; x ≤ Γx (t) T (x,t) = C1 The mathematical model consists of (15) – (18) in dependent β erf variables e, q, and T , or (17) – (19) in dependent variables e and T . 2 (24) In the published works specific heat c p is generally considered as a     β x function of temperature, but in general ρ = ρ(T ), c p = c p (T ), and erf − erf √ 2 2 t + t0 k = k(T ) are permissible but can only be used outside the transition T (x,t) = C   ; x > Γx (t) 2 β region. Consider equation (17) during phase change, i.e. change in erfc s 2 e from es to el . When α = e−e L f and T = Ts , (17) is identically satisfied. α = 1 for e > es + L f clearly indicates instantaneous addition The interface location Γx (t) is defined by of latent heat. Both models in e, q, T and e, T have been used in √ the published works [10–14]. Γx (t) = β t + t0 (25) The parameter β is obtained by solving the equation  

Model (b) If we assume that c p , k, and ρ are constant in the solid and liquid regions and have values c ps , ks , ρs and c pl , kl , ρl , then we can write explicit forms of (19) for solid and liquid phases by using es = c ps T and el = c pl T . These equations are augmented by a heat 68

The phase field mathematical models of phase change also introduce a finite width variable transition region between the two states. These models are based on the work of Cahn and Hilliard [4] and are derived using Landau-Ginzberg theory of critical phenomena [5]. A phase field variable p is introduced which has a value of −1 for solid phase and +1 in the liquid phase. The length of (2) The sharp interface creates singularity of e at T = Ts which the transition region between the solid and the liquid phases is conposes many obvious difficulties in the computation of the nu- trolled by choosing a value of ξ (figure 2) that corresponds to inmerical solutions of the associated initial value problem. termediate value of p. (1) One of the major disadvantages of the sharp-interface mathematical models is that the phase change is assumed to occur at a constant temperature Ts . Thus, e changes from es to el at constant T = Ts . This is true regardless of the form of the mathematical models.

(3) It is meritorious to eliminate q as a dependent variable as done in case of (19) as it reduces the number of dependent variables in the mathematical model. But this reduction is at the cost of appearance of the second derivative of T with respect to spatial coordinates in the energy equation, which in context of finite element methods of approximation requires higher order regularity for the approximation of T .

(4) In addition to eliminating q as a dependent variable, the specific internal energy e can also be substituted in the energy equation yielding a single nonlinear diffusion equation in temperature T . We postpone details of this until a later section. (5) It is critical to point out that all sharp-interface models are derived based on a priori existence of the transition front as initial condition. As a result these models cannot simulate initiation of the transition front. The models simply simulate propagation of this front during evolution. This is vital physics that is necessary in almost every phase change application and is missing in the sharp-interface approach.

Figure 2: Expected spatial profile of phase field through a solid-liquid interface

The phase field approach avoids the explicit treatment of the interface conditions as employed in the sharp-interface models. Instead we use a coupled system of nonlinear evolution equations in temperature T and phase variable p [24]

(6) Finally, if one considers computations of the numerical solutions for phase change processes to be essential, then sharpinterface model of phase change processes are not meritorious.

In smooth-interface models the phase change is assumed to take place over a finite temperature range [Ts , Tl ] (see figure 1(b)) during which e is continuous and differentiable in temperature T. The range [Ts , Tl ], referred to as transition region consisting of solidliquid mixture i.e. a mushy region, can be as narrow or as wide as desired. At T = Ts the state of the matter is solid whereas at T = Tl it is pure liquid. Since the properties ρ, c p , k have different values for solid and liquid phases, it is often meritorious to consider these as functions of temperature T with continuous and differentiable behavior for Ts ≤ T ≤ Tl between their values ρs , c ps ,ks and ρl , c pl , kl for solid and liquid states respectively. In the following we present details of two smooth-interface mathematical models, one based on phase field approach and the other based on the energy equation (12) with transition region [Ts , Tl ] in which ρ, c p , k and L f are continuous and differentiable functions of the temperature T .

in which α is related to the kinetic parameter [24], ∆ = ∂∂x2 + ∂∂y2 + ∂2 , and f = f (p, T ) is referred to as the restoring potential or free ∂ z2

energy potential. Equations (27) and (28) can be interpreted in a simple way. Equation (28) is a linear time evolution of p governed by imbalance between the excess interface free energy and the restoring potential f (p, T ). The energy equation (27) has a source term 12 L f ∂∂tp to account for the latent heat release or absorption at the moving interface. When the phase field equations (27) and (28) are employed to simulate real solidification or melting problems, we expect that sharp-interface conditions are approached as the interface thickness ξ → 0. The results in phase field models unfortunately depend largely on thermodynamic consistency of the potential f (p, T ). The work of Caginalp [25] provides a strong indication that the sharp-interface limit is attained for all forms of

free energy potential f (p, T ) in which T and p coupling is lin∂ 2 f (p,T ) ∂ p∂ T = 0.

ear i.e. Given a specific form of f (p, T ), the entropy/energy/temperature scales must obey the relationships: ∆η = (ηliquid − ηsolid ) =

(2) As in case of sharp-interface models, here also the interface must be defined as initial condition. The phase field models are not capable of initiating phase transition. Obviously, this is a major drawback of these models. This drawback is due to the use of double well function f (p, T ).

(3) When the spatial domain is either solid or liquid, the free energy density functions used presently do not allow initiation of the transition zone or front due to the presence of two distinct in which η is entropy, ∆η is change in entropy, and Tm is mean temminima, regardless of the temperature. For example if the spaperature. In the Caginalp Potential (CP) model [24], dependence tial domain is liquid and heat is removed from some boundary, of f on T is taken into account by adding a simple linear term to the liquid will remain in the liquid state although the temperathe double well potential in p. ture may have fallen below the freezing temperature. ∆η 1 2 2 pT (31) f (p, T ) = (p − 1) − 8a 2 (4) This drawback of phase field models presents serious problems The parameter a is chosen such that ∂∂ pf exhibits three distinct in simulating phase transition processes in which initiation and detection of the location of the transition zone is essential as it roots, near 0 and ±1. From (31) we note that minima of f (p, T ) at may not be known a priori. p = ±1 changes as T departs from zero. Figure 3 shows a plot of p versus f (p, T ) for T = 0, T < 0, and T > 0 (with a = ∆η = 1). ∆η T =0 =

(5) When the phase transition region is specified as initial condition, the phase field models predict accurate evolution i.e. movement of the transition region.

Mathematical models used in the present work: smoothinterface models (first group of models)

The mathematical models used in the present work presented in this section are derived based on the assumptions that the transition region between the liquid and solid phases occurs over a small temperature change (width of the transition region [Ts , Tl ]) in which specific heat, thermal conductivity, density, and latent heat of fusion and hence specific internal energy change in a continuous -1.5 -1 -0.5 0 0.5 1 1.5 Phase Variable, p and differentiable manner. Figures 4(a),(b),(c),(d),(e) show distributions of ρ, c p , k, L f , and e in the transition region [Ts , Tl ] between Figure 3: Double-well behavior of restoring potential for various values the solid and liquid phases. The range [Ts , Tl ] i.e. the width of the of temperature transition region, can be as narrow or as wide as desired by the physics of phase change in a specific application. The transition region is assumed to be homogeneous and isotropic. This assumpFor a finite value of a, a small amount of latent heat is retion is not so detrimental as in this case the constitutive theory only leased at positions away from the interface. This undesirable efconsists of the heat vector due to the zero velocity field and zero fect fades as a → 0. Indeed, Caginalp et al. [25–28] have estabstress assumptions. lished that as ξ → 0 and a → 0 the phase field equations (27) and The mathematical models derived and presented here are same (28) with f (p, T ) defined by (31) produce solutions that approach in Lagrangian as well as Eulerian description, and are based on the sharp-interface limits. first law of thermodynamics using specific total energy and the heat vector augmented by the constitutive equation for the heat vector Remarks (Fourier heat conduction law) and the statement of specific total (1) The phase field models require free energy potential f (p, T ). energy incorporating the physics of phase transition in the smooth There are some guidelines to establish this but for the most part interface zone between liquid and solid phases. the procedure is not deterministic. 70

If Q(T ) represents any one of the quantities ρ(T ), c p (T ), k(T ), and L f (T ), then we define  ; T < Ts  Qs Q(T ) ; Ts ≤ T ≤ Tl Q(T ) = (37)  Ql ; T > Tl

when n = 3, Q(T ) is a cubic polynomial in T . The coefficients c0 and ci , i = 1, 2, 3 in (38) are calculated using the conditions:

when n = 5, Q(T ) is a 5th degree polynomial in T . The coefficients c0 and ci , i = 1, ..., 5 in (38) are calculated using the conditions:

Figure 4: ρ, c p , k, L f and e in the smooth interface transition region between the solid and liquid phases as functions of Temperature T

Remarks (1) By letting Q to be ρ, c p , k and L f , dependence of these properties on temperature can be easily established.

In the absence of sources and sinks we have De ∇·q = 0 +∇ Dt

(32) (2) In case of L f we note that L f (Ts ) = 0 and L f (Tl ) = L f (value of latent heat of fusion). Assuming Fourier heat conduction law as constitutive theory (3) Thus all transport properties including latent heat of fusion are for q , we can write explicitly defined as functions of temperature T in the transition region. ∇T q = −k(T )∇ ∀(xx ,t) ∈ Ωx × Ωt = Ωx × (0, τ) (33) ρ

 ∂T ∂T ∂Lf c p (T )dT + L f (T ) = c p (T ) + (35) ∂t ∂t ∂t T0 and

∂ L f (T ) ∂T ∇·q = 0 +ρ(T ) +∇ ∂t ∂t ∀(xx ,t) ∈ Ωx × Ωt = Ωx × (0, τ)

Hence, (36) can be written as   ∂ L f (T ) ∂ T ∇·q = 0 ρ(T )c p (T ) + ρ(T ) +∇ ∂T ∂t

Equations (42) and (43) are smooth-interface mathematical (36) model in dependent variables T and q. ρ(T ), c p (T ), k(T ), and L f (T ) are defined using (37)-(40). By substituting q from (43) into 71

dependent variable and use L f (T ) = G(T ) as addition equation in which G(T ) is functional relationship of L f on T .

(42), we obtain a single nonlinear diffusion equation for smooth interface phase change model.  
∂ L f (T ) ∂ T ∇· k(T )∇ ∇ (T ) = 0 (44) −∇ ρ(T )c p (T ) + ρ(T ) ∂T ∂t

  ∂Lf ∂T ρ(T )c p (T ) + ρ(T ) ∂ T ∂t   ∂ k(T ) 3 ∂ T 2 − ∑ ∂ xi − k(T )∆T = 0 ∂ T i=1

or   ∂ L f (T ) ∂ T ρ(T )c p (T ) + ρ(T ) ∂T ∂t 2 3  ∂ k(T ) ∂T − ∑ ∂ xi − k(T )∆T = 0 ∂ T i=1 2

(45) This model is second order in T but first order in L f . Model D: Dependent Variables T , q , and L f In this case we consider the mathematical model (47), but also introduce L f as a dependent variable.

where ∆ = ∂∂x2 + ∂∂x2 + ∂∂x2 . x1 = x, x2 = y, and x3 = z have been used for convenience. Equation (45) is the final form of the mathematical model in temperature T .

(1) The mathematical models presented in this section can be written in alternate forms. These are summarized in the following based on choice of dependent variables.

This is a first order model in T , q, and L f . Model A: Dependent Variable T If we consider T as the only dependent variable, then the mathematical model is given by (45) i.e.     ∂ L f (T ) ∂ T ∂ k(T ) 3 ∂ T 2 − ρ(T )c p (T ) + ρ(T ) ∑ ∂ xi (46) ∂T ∂t ∂ T i=1 − k(T )∆T = 0

(2) The mathematical models given in remark (1) are all valid models. We present more discussion on these models in the section on numerical studies.

Second group of mathematical models for phase change based on nonzero stress and velocity fields in all phases

This model requires higher order regularity of approximations of T in finite element processes of calculating numerical soWhen the media are not stress free and the velocity field is lutions for T . This is due to second order derivatives of the not zero, the mathematical models for phase change processes retemperature with respect to spatial coordinates appearing in quire use of all conservation and balance laws for solid and liquid (46). phases, as well as the transition region. The mathematical models must incorporate the physics of solid, liquid, and transition regions Model B: Dependent Variables T , q and their interactions during the evolution of the phase change proIn this case the mathematical model consists of equations (42) cess. In the approach discussed here the mathematical models for and (43). all phases are strictly based on conservation and balance laws and    the transition region is assumed to be a smooth interface between ∂Lf ∂T   the solid and the liquid phases. We consider details of the models ∇·q = 0 ρ(T )c p (T ) + ρ(T ) +∇ ∂ T ∂t (47) for all three phases and present discussion regarding their validity   q = −k(T )·· ∇ T and use in determining phase change evolution. This is a system of first order partial differential equations in T and q, hence lower order regularity on both q and T compared 2.2.1 Liquid Phase to T in Model A. If we assume (for simplicity) the liquid phase to be incompressible Newtonian fluid with constant properties, then the mathematModel C: Dependent Variables T , L f ical model for this phase is standard continuity, momentum equaIn the mathematical model, rather than replacing L f (T ) with tions, energy equation, and the constitutive theories for contravarian expression, a function of T , we could also consider L f as a ant deviatoric Cauchy stress tensor and heat vector in Eulerian de72

scription with transport. In the absence of body forces, we have  ρ̄l ∇¯ · v̄v = 0        (0)    ∂ v̄i ∂ v̄i ∂ p̄ ∂d σ̄i j   ρ̄l + v̄ j + − =0   ∂t ∂ x̄ j ∂ x̄i ∂ x̄ j     (50) ∂ T̄ (0) ¯ ¯  + v̄v · ∇ T̄ + ∇ · q̄q − d σ̄ ji D̄i j = 0 ρ̄l c̄ pl   ∂t     (0)    d σ̄i j = 2 µ̄ D̄i j     ¯ q̄q = −k̄l ∇ T̄

mechanical pressure p to d σ . With this decomposition this mathematical model appears to have the same dependent variables (but not necessarily the same physical meaning) as the one for fluid in Section 2.2.1. Hypoelastic Solid

p̄ is mechanical pressure assumed positive when compressive. ρ̄l , c̄ pl , k̄l , µ̄ are the usual constant transport properties of the medium. We remark that x̄x are fixed locations at which the state of the matter is monitored as time elapses i.e. x̄x location is occupied by different material particles for different values of time. In this mathematical model material point displacements are not monitored. 2.2.2

In the solid phase the most appropriate form of the mathematical model can be derived using conservation and balance laws in Lagrangian description. Hyperelastic Solid If we assume the solid phase to be hyperelastic solid matter, homogeneous, isotropic, and incompressible with infinitesimal deformation and constant material coefficients, then we have the following for continuity, momentum equations in the absence of body forces, energy equation, and the constitutive equations (using σ for stress tensor).  ρ0 = ρs as |J| = 1      ∂ vi σi j    − =0 ρs   ∂t ∂xj      ∂T   ∇ · q ρs c ps +∇ = 0    ∂t  σi j = Di jkl εkl ∀(xx ,t) ∈ Ωx t = Ωx × Ωt (51)       1 ∂ ui ∂ u j   εi j = +   2 ∂ x j ∂ xi      ∂ ui    vi =   ∂t    q = −ks∇ T

If we assume the solid phase to be hypo-thermoelastic solid matter, isotropic, homogeneous, and incompressible with constant material coefficients then the mathematical model can be derived in Eulerian description with transport. The constitutive theory for the stress tensor for such materials is a rate theory of order one in stress and strain rate tensors i.e. convected time derivative of order one of the stress tensor is related to the convected time derivative of order one of the conjugate strain tensor. If we consider σ (0) decomposition then we have the following for σ (0) = − p̄II + d σ̄ σ̄ continuity, momentum and energy equations, and the constitutive equations.  ρ̄s∇¯ · v̄v = 0        (0)    ∂ v̄i ∂ p̄ ∂d σ̄i j ∂ v̄i   + v̄ j − =0 + ρ̄s   ∂t ∂ x̄ j ∂ x̄i ∂ x̄ j     (52) ∂ T̄ ¯ T̄ + ∇ ¯ · q̄q = 0  + v̄v · ∇ ρ̄s c̄ ps   ∂t     (1) (1)    d σ̄i j = D̃i jkl γkl     ¯ q ∇ q̄ = −k̄s T̄ σ (1) is the first convected time derivative of the deviatoric conσ̄ travariant Cauchy stress tensor and γ (1) is the first convected time derivative of the Almansi strain tensor, a contravariant measure of strain. It has been shown [29] that for thermo-hypoelastic solids the σ (0) ) = 0 in R3 , continuity equation in (52) must be replaced by tr(d σ̄ (0) (0) σ ) − p̄ = 0 in R2 , and tr(d σ̄ σ ) − 2 p̄ = 0 in R1 as additional tr(d σ̄ σ (0) . We note that in equation relating mechanical pressure p̄ to d σ̄ hypoelastic solids, strain rate produces stress as opposed to strain as in the case of hyperelastic solids. Secondly, such a model allows transport which is not present in the deformation of thermoelastic solids. 2.2.3 Transition Region In the transition region the consideration of the physics of phase transformation and how we account for it in the development of the mathematical model determines the ultimate outcome of the details of the mathematical model. The following approaches are used or are possibilities.

In this description the locations x are locations of material points, hence the deformation of the material points is monitored during evolution. We note that we can also introduce the stress decomposition σ = −pII + d σ with tr(d σ ) = 0 in R3 , tr(d σ ) − p = 0 in R2 , and tr(d σ ) − 2p = 0 in R1 as additional equation relating 73

(a) We can assume the transition region as a homogeneous, saturated mixture of fluid and solid constituents with appropriate volume fractions based on temperature. In this approach the solid particles are always mobile, which poses problems as we approach the solid phase. The choice of Lagrangian or Eulerian description (with transport) is also not straightforward.

This approach has not been used in phase transition applications.

(b) We assume that freezing or melting in the transition region creates a porous media with variable permeability. This approach has been used but not in conjunction with the full NavierStokes equations.

Third group of mathematical models for phase change: the stress field is assumed to be constant and velocity field is assumed zero in the solid phase but both are nonzero in the liquid phase and transition region

The mathematical models in section 2.2 fail to provide interaction between the phases. We note that the main source of this (c) Some variations of mixture theory with various approximaproblem is that not all dependent variables in the two descriptions tions are possible. describe the same physics. For example velocities in Lagrangian description are time rate of change of displacements of a material point, whereas in Eulerian description the velocities at a location Remarks are velocities of different material points for different values of It is perhaps more straightforward to illustrate the problems as- time. Since we want to consider phase transition in the presence sociated with these mathematical models and their use in phase of flow, the mathematical models for fluid derived based on contransition if we consider sharp interface between solid and liquid servation and balance laws must remain intact. In the solid region regions. In this case purely solid phase is in contact with purely we have displacements, their time derivatives (velocities) and stress (dependent on displacements) in Lagrangian description. These liquid phase. We consider the following: quantities do not have the same physical meaning in case of fluid (1) Lagrangian description with hyperelastic solid assumption is using Eulerian description with transport, thus must be eliminated ideal for the solid phase and the Eulerian description with if we seek interaction the two media without changing the mathetransport is suitable for the liquid phase. In the Lagrangian matical model for fluid. This gives rise to constant stress field and description the locations x are the positions of the material par- zero velocity field in the solid region. In the transition zone the veticles that undergo evolution and thus we have displacements locity field must transition from nonzero state at the liquid boundary to zero state at the solid boundary, and the stress must assume of each material particle in time during evolution. a constant value. Thus in this approach only the energy equation On the other hand, in Eulerian description with transport the provides the connecting link between the solid and transition relocations x̄x are fixed locations that are occupied by different gions. In the solid region, the energy equation has no transport material particles for different values of time. Thus, in this ap- terms (due to Lagrangian description) and no dissipation terms as proach we do not have displacement history of each material the solid phase is thermoelastic, but these would have been zero particle in time during evolution. even otherwise as the velocity field and divergence of the stress field are zero. In the liquid region we have energy equation with At the interface between the solid and the liquid regions, these transport as well as dissipation, both of which approach zero in the two mathematical models do not provide interaction. This has transition region as the state evolves from liquid to solid and hence been established by Surana et al. [30]. Forcing these matheyields the desired energy equation for the solid phase. matical models to interact will produce spurious behavior. Thus, in this approach we assume that the solid phase has constant stress field and the velocity field is zero in this phase but in (2) We could consider hypoelastic solid description (Eulerian dethe liquid phase we consider full Navier-Stokes equations based on scription with transport) for the solid phase and the Eulerian conservation and balance laws. Some aspects of the approach disdescription with transport for the liquid phase. In this case, cussed here are also found in [17, 19–21] but differ significantly in interaction between the two phases is intrinsic in the mathethe specific details of the mathematical model and numerical commatical model and is mathematically consistent. putations of the evolution. We consider that the solid and liquid However the hypoelastic solids have transport which is not the phases have smooth interface in which all transport properties vary case for solid phase and the first order rate constitutive theory in a continuous and differentiable manner as in Section 2.1.2. We is nonphysical for solid phase as it could yield zero stress in consider details of the mathematical models in the following. the absence of strain rates. The presence of transport for the solid phase is also problematic during the liquid to solid phase change. Thus when the stress field and the velocity field are not zero the current mathematical models for solid and liquid phases do not permit interaction of the solid and liquid phases (see Surana et al. [30]).

Liquid Phase For this phase we consider standard Navier-Stokes equations in Eulerian description (with transport) as used for fluids (constitutive 74

theory for stress based on Newton’s law of viscosity). ρ̄l ∇¯ · v̄v = 0   (0) ∂ v̄i ∂ v̄i ∂ p̄ ∂d σ̄i j + v̄ j ρ̄l + − =0 ∂t ∂ x̄ j ∂ x̄i ∂ x̄ j   ∂ T̄ (0) ρ̄l c̄ pl + v̄v · ∇¯ T̄ + ∇¯ · q̄q − d σ̄ ji D̄i j = 0 ∂t (0) d σ̄i j = 2 µ̄ D̄i j

The mathematical models for solid, liquid, and transition phases can be combined into a single mathematical model. f¯l ρ̄(T̄ )∇¯ · v̄v = 0   (0) ∂ σ̄ ¯fl ρ̄(T̄ ) ∂ v̄i + v̄ j ∂ v̄i + f¯l ∂ p̄ − d i j = 0 ∂t ∂ x̄ j ∂ x̄i ∂ x̄ j    ∂ L̄ f (T̄ ) ∂ T̄ + f¯l v̄v · ∇¯ T̄ ρ̄(T̄ ) c̄ p (T̄ ) + ∂t ∂ T̄
(0) + ∇¯ · q̄q − f¯l d σ̄ ji D̄i j = 0 (0)
f¯l d σ̄i j = 2µ̄ D̄i j

Solid Phase Since in the solid phase the stress field is assumed constant and the velocity field is assumed zero, the mathematical model for this phase only consists of the energy equation and the constitutive theory for heat vector. In the absence of the velocity field and stress field, there is no distinction between the Lagrangian and the Eulerian descriptions, but we use overbar to provide transparency between this description and the one given by (53). σ (0) = 0 ∇¯ · d σ̄ ∂ T̄ + ∇¯ · q̄q = 0 ρ̄s c̄ ps ∂t q̄q = −k̄s∇¯ T̄

From (55) and (56) we note that in the solid phase f¯l = 0, hence ∂ L̄ f ∀(x̄x ,t) ∈ Ωx̄x t = Ωx̄x × Ωt (54) continuity equation is identically zero, ∂ T̄ = 0 and the others reduce to (0)

In this region we note that the momentum equations in (53) must be satisfied for zero velocity field with zero pressure gradient. From the first set of equations in (54) we note that a constant deviatoric Cauchy stress field is admissible. The values of the constant stresses in the solid region are determined by the values of the stresses at the solid-liquid interface. In other words, the constant of integration in the first set of equations in (54) is determined using the values of the stresses at the solid-liquid interface. Since the stress gradients are zero in the solid region, the stress values in the solid region remain constant and their values are same as those at the liquid-solid interface (see numerical studies in section 3.5).

Comparing (57) with (53) we note that v̄ = 0 and thus D̄i j = 0, ∂d σ̄

and ∂ x̄ijj = 0, therefore (53) reduces to (57), which is same as (54). Presence of the first equation in (57) is essential as it ensures ∂d σ̄

that it is oscillation free so that ∂ x̄ijj = 0 would hold precisely Transition Region everywhere in the solid phase. Thus the challenge in the matheIn the transition region from liquid to solid the mathematical matical model (55),(56) is to ensure that D̄i j = 0 is achieved in the ensure that v̄ and its gradients as well as model transitions from (53) to (54) or vice versa. As in Sec- solid phase which would (0) σ the gradients of σ̄ are zero in the solid phase. The momend tion 2.1.2, we consider a transition region [T̄s , T̄l ] in temperature. tum equation in (55) (or (57)) when satisfied for zero velocity field In this region we assume that k̄, c̄ p , ρ̄ transition from solid to liquid (0) ∇ · σ ensures that σ̄ is identically zero and oscillation free in the d values in a continuous and differentiable manner as described in solid phase. ¯ ¯ Section 2.1.2. Let fl and fs be the liquid and solid fractions with f¯s = 1 − f¯l and 0 ≤ f¯l ≤ 1 in the transition region. We also asThis mathematical model is used in the present work to present sume that release or absorption of latent heat of fusion L f is also numerical studies for phase change when the stress field and the continuous and differentiable in the transition region. velocity field in the liquid phase are not zero. 75

3.1. Dimensionless form of the mathematical

models used in the present work In the following we present dimensionless form of the mathematical models of phase change based on: (i) the assumption that stress field and velocity field are zero in solid, liquid, and transition phases, and (ii) the assumption that in the solid phase the stress field and velocity field are zero but in the liquid phase full NavierStokes equations constitute the mathematical model. In both models the transition zone of width [Ts , Tl ] in temperature is assumed homogeneous and isotropic in which ρ, c p , k make transition from solid to liquid phase and vice versa in a continuous and differentiable manner. In order to nondimensionalize the mathematical models we choose reference quantities to obtain dimensionless dependent and independent variables and other quantities. The quantities with hat ( ˆ ) are with their usual dimensions, quantities with zero subscript are reference quantities and the quantities without hat ( ˆ ) are dimensionless quantities. We define xi = x̂i /L0

Since the velocity field is assumed zero, t0 cannot be defined using L0 and v0 . We can choose the following: t0 = L02 ρ0 c p0 /k0

Using (66), the mathematical model (64) and (65) reduces to ∂T ∇· q +ρ ρc p +∇ ∂t ∇T q = −k∇

Equations (67) and (68) are a system of first order PDEs in T and q in which reference time t0 and reference latent heat of fusion L f 0 are defined by (66). Alternatively, if we substitute q from (68) into (67), then we obtain a single PDE in temperature T .

    Lf0 ∂Lf t0 k 0 ∂T ∇ · q + + ρ =0 2 ∂t c p0 T0 ∂t L0 ρ0 c p0

Equation (69) contains up to second order derivatives of temperature T in space coordinates. The mathematical models (67) and 3.1.1 Mathematical model based on the assumption of (68) as well as (69) can be used in numerical studies, but the choice zero stress and zero velocity field in all phases (first of local approximations for minimally conforming approximation group of models) ∂L spaces differ in the two. Since L f = L f (T ), ∂ Tf is strictly deterRecall the following mathematical model presented in Section ministic. Other mathematical models (Model C and Model D) presented in Section 2.1.2 have similar dimensionless forms. 2.1.2. 76

(0) Mathematical model when the stress field is asfield i.e. in the solid phase ∂ v̄i /∂ x̄ j = 0 and ∂ d σ̄i j /∂ x̄ j = 0 sumed constant and the velocity field is assumed must hold in the solid phase. zero in the solid phase but nonzero in both the liq(2) Based on (1), it may be possible to redefine new dependent uid and transition regions (third group of models) variables so that during numerical computations, conditions in Recall the mathematical model given by (55) and (56) (1) are also satisfied with this choice. This indeed is the case as shown in the model problems in section 3.5.1.   f¯l ρ̄ˆ (T̄ˆ )∇ˆ¯ · v̄vˆ = 0        (0)  ˆ

3.2. Computational methodology for computing

∂ σ̄ ˆ ˆ ˆ  ∂ v̄ ∂ v̄ ∂ p̄ d i j i i   f¯l ρ̄ˆ (T̄ˆ ) + v̄ˆ j + f¯l − =0  evolution of IVP describing phase change  ∂ tˆ  ∂ x̄ˆ j ∂ x̄ˆi ∂ x̄ˆ j     ˆ     The mathematical models describing phase change are a system ∂ L̄ˆ f (T̄ˆ ) ∂ T̄ + f¯l v̄vˆ · ∇ˆ¯ T̄ˆ ρ̄ˆ (T̄ˆ ) c̄ˆ p (T̄ˆ ) + (70) of nonlinear partial differential equations. Numerical solutions are ∂ tˆ  ∂ T̄ˆ   computed using space-time least squares finite element processes
 (0)  + ∇ˆ¯ · q̄qˆ − f¯l d σ̄ˆ ji D̄ˆ i j = 0   for a space-time strip (in R1 ) or a space-time slab (in R2 ) the with   
 time marching. The mathematical models utilized in the computa(0)   f¯l d σ̄ˆ i j = 2µ̄ˆ D̄ˆ i j   tional studies are a system of PDEs. In case of R1 , the space-time    domain of a space-time strip for an increment of time is discretized q̄qˆ = −k̄ˆ (T̄ˆ )∇ˆ¯ T̄ˆ using nine-node p-version space-time elements. In case of R2 , the where space-time slab is discretized using 27-node p-version space-time f¯l = 1 ; liquid phase elements. Local approximations of class C0 and C1 in space and f¯l = 0 ; solid phase (71) time are in the computations. 0 ≤ f¯l ≤ 1 ; T̄ˆs ≤ T̄ˆ ≤ T̄ˆl For an increment of time i.e. for a space-time strip or a slab, solution of the non-linear algebraic systems is obtained using NewDimensionless forms of (70) and (71) can be obtained using ton’s linear method with line search. Newton’s linear method is (58): considered converged when the absolute value of each component  f¯l ρ̄(T̄ )∇¯ · v̄v = 0 of δ I = {g} is below a preset threshold ∆, numerically computed    −6 has been used in all numerical studies. Discretiza     zero. ∆ ≤ 10   ∂ v̄i p0 ∂ p̄ ∂ v̄i   + v̄ j + f¯l f¯l ρ̄(T̄ ) tion and p-levels (considered to be uniform in space and time) are   ∂t ∂ x̄ j  ρ0 v20 ∂ x̄i  chosen such that the least squares functional I resulting from the     (0)   ∂d σ̄i j residuals for the entire space-time strip or slab is always of order  τ0  =0  −  2 of O(10−6 ) or lower and hence good accuracy of the evolution is  ∂ x̄ ρ0 v0  j   always ensured.     L f 0 ∂ L̄ f (T̄ ) 1 ∂ T̄ ρ̄(T̄ ) c̄ p (T̄ ) + 2 + f¯l v̄v · ∇¯ T̄  (72) Ec ∂t  v0 ∂ T̄    

3.3. 1D phase change model problems

 
1 ¯ τ  (0) 0   + ∇ · q̄q − f¯l σ̄ D̄ = 0 ij d ji  We consider three model problems. In the first model problem  ReBr ρ0 v20      we present a comparison of the smooth-interface solutions (present  
  ¯fl d σ̄ (0) = µ0 v0 2µ̄ D̄i j approach) with the theoretical solution obtained using the sharp  ij  L0 τ0  interface method. In the other two model problems we consider    ¯ solid-liquid and liquid-solid phase change. q̄q = −k̄(T̄ )∇ T̄

f¯l = 1 f¯l = 0 0 ≤ f¯l ≤ 1 ρ0 v0 L0 Re = µ0 µ0 v20 Br = k0 T0

Model Problem 1: Comparison of Sharp- and Smooth-Interface Solutions

The sharp interface solution [12] has only been reported for (73) constant material coefficients. When the material coefficients vary, ; Reynolds Number i.e. are a function of temperature, the theoretical solution of the resulting mathematical model has not been reported, perhaps due ; Brinkman Number to complexity. We choose ρ = 1, c p = 1, k = 1, and L f = 1. The spatial domain consists of 0 ≤ x ≤ 1. Figure 5 shows a space-time strip Ωxt = [0, 1]×[0, ∆t]. The space-time domain Ωxt is discretized Remarks using a uniform mesh of 500 p-version nine node space-time ele(1) We keep in mind that in the solid phase the momentum equa- ments. The spatial domain [1, 4] × [0, ∆t] is discretized using a 30 tions must be satisfied for zero velocity field and constant stress element uniform mesh. The spatial domain 1 ≤ x ≤ 4 is added to 77

0 ≤ x ≤ 1 to approximate the boundary condition at x = ∞ in the theoretical solution with x = 4 in the computed solution. t

In the smooth-interface solutions we also use ρ = 1, c p = 1, k = 1, and L f = 1 i.e. constant material coefficients regardless of phase. The evolution is computed using ∆t = 0.01 for 100 time steps i.e. up to t = 1.0. The latent heat L f is expressed as a polynomial in temperature T in the transition zone. Generally a cubic or fifth degree polynomial in T for L f is found adequate (equations (37)–(40)). Evolution of temperature and latent heat for 0 ≤ t ≤ 1 from smooth interface and comparison with sharp-interface solution are shown in figures 7 and 8. Interface location versus time t from smooth and sharp interface locations are compared in figure 9. Center of the transition region is considered as interface location in smooth interface solution.

Figure 5: Schematic of first space-time strip, BCs, IC, and spatial discretization

The initial conditions on temperature T at t = 0 are defined piecewise by the following. ;

√
erf β /2 − erf x/2 t0
Θ(x) = C1 erf β /2
√
erf β /2 − erf x/2 t0
Θ(x) = C2 erfc β /2

From figures 7–9 we note that smooth-interface solutions are in good agreement with sharp-interface solutions. The sharp-interface theoretical solution is only possible for constant ρ, c p , and k, whereas smooth-interface solutions are possible for variable ρ, c p , and k. Smooth-interface solutions with transitions in material coefficients due to phase change describes physics of phase transitions more precisely.

In the theoretical solution for sharp-interface (23)–(26), the following coefficients are used. C1 = −0.085 t0 = 0.1246

The mathematical model (69) is used for computing smooth0 0.2 0.4 0.6 0.8 1 interface solutions. For smooth-interface solutions the transition Distance, x region is defined by [Ts , Tl ] = [−0.001, 0.001]. p-levels in space Figure 7: Model Problem 1: Evolution of temperature using smoothand time are chosen to be 7, with solutions of class C1 in space and interface model and sharp-interface theoretical solution, time. Figure 6 shows a plot of the initial condition at t = 0. C11 (Ω̄ext ), p = 7, ∆t = 0.01 0.02 1.2

Figure 6: Initial condition Θ(x) at t = 0, temperature distribution from the theoretical solution of the sharp-interface model

Figure 8: Model Problem 1: Evolution of latent heat (smooth interface), C11 (Ω̄ext ), p = 7, ∆t = 0.01

Model Problem 2: 1D Liquid-Solid Phase Change; Initiation and Propagation of Phase Transition

In this model problem we consider 1D liquid-solid phase change with variable material coefficients and to demonstrate the Figure 9: Model Problem 1: Interface location as a function of time, ability of the proposed formulation in initiating phase transition as C11 (Ω̄ext ), p = 7, ∆t = 0.01 well as in simulating its evolution as time elapses. Numerical solutions are calculated and compared for constant density (ρ̂ = ρ̂l in all phases) as well as variable density. Figure 10 shows space-time 3.3.2 Transport properties and reference quantities for strip Ω̄xt = [0, 1] × [0, ∆t], initial conditions, and boundary condiliquid-solid and solid-liquid transition numerical tions. studies with zero stress and velocity fields in all tˆ phases 0

In all numerical studies using zero velocity and zero stress field for the entire domain, we consider the liquid phase to be water and the solid phase to be ice with the following properties.

In the transition region ρ(T ), c p (T ), k(T ) and L f (T ) are ast = ∆t sumed to vary in a continuous and differentiable manner between (c) Boundary condition dT dx at x = 1 the temperatures Ts and Tl defining the transition region between solid and liquid phases. Figure 10: Liquid-solid phase transition: space-time strip, boundary conditions, and initial condition

Reference quantities: Regardless of solid-liquid or liquid-solid phase transition we consider the following reference quantities:

We consider solutions of class C11 with p-level of 9 in space and time. With this choice the space-time integrals are Riemann 79

Figures 12(a), 14(a), and 15(a) show differences in the evolution of L f , c p , and k for constant and variable densities even in the very early stages of the evolution. Constant density results lag variable density solutions.

in time but Lebesgue in space. This choice functions quite well in simulating the evolution (low residuals). We choose the phase transition zone [Ts , Tl ] to be [−0.001, 0.001]. A different (smaller or larger) choice of transition zone width in temperature does not alter the location of the center of the transition region.

Figure 11: Model Problem 2: Evolution of temperature for liquid-solid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Computed numerical results are presented in figures 11–15. Figures 11(a), 12(a), 13(a), 14(a), 15(a) show plots of T , L f , ρ, c p , and k versus x during initial stages of the evolution (0 ≤ t ≤ 0.2). Continuous extraction of heat from the right boundary progressively lowers the temperature at the boundary and in the neighborhood of the boundary which eventually results in the initiation of phase change. Variations in L f (T ), c p (T ), k(T ) and ρ(T ) follow changes in temperature during evolution. From figure 11(a) we note that both constant and variable densities yield almost the same evolution of the temperature during initial stages of the evolution.

Figure 12: Model Problem 2: Evolution of latent heat for liquid-solid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figures 11(b), 12(b), 13(b), 14(b), 15(b) show fully formed phase change transition region (liquid to solid) beginning with t = 0.8 and its propagation during evolution (0.8 ≤ t ≤ 4.8). For most space-time strips during time marching using ∆t = 0.04, I < O(10−6 ) and |(gi )|max ≤ 10−6 ensure accurate evolution that satisfies GDE quite well over the entire space-time domain of each space-time strip. Evolutions of all quantities are smooth and free of oscillations. The influence of variable density can be seen clearly in these graphs. The variable density results lead constant density evolution and the difference between them increases as the evolution proceeds. 80

From figure 11(b) we clearly observe linear heat conduction in The differences in the computed solutions for constant and variliquid and solid phases (constant but different slopes of T versus x) able density are noticeable. We note that the center of the phase separated by smooth transition region. transition zone for variable density case is ahead of the constant density case during the entire evolution and the difference between the two increases as the evolution proceeds. 1.1

Figure 13: Model Problem 2: Evolution of density for liquid-solid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figure 14: Model Problem 2: Evolution of specific heat for liquid-solid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

If we define the center of the transition zone as the location x of the phase front, then using the results in figures 11–15 we can plot a graph of location x versus time t marking the location of the phase change front in time. Figure 16 shows such a plot for the results presented in figures 11–15. The transition region width for these numerical studies consist of [Ts , Tl ] = [−0.001, 0.001]. It is also obvious from figures 11–16 that the choice of constant density in all phases (ρ = ρl used here), as is commonly used in the published works, will produce results that do not agree with the actual physics of phase change (variable density).

Similar studies were repeated for [Ts , Tl ] = [−0.002, 0.002] i.e. double the width of the transition zone, with virtually no change in the location of the center of the transition region. The phase transition evolution for this model problem cannot be simulated using sharp-interface and phase field approaches as this model problem requires initiation of phase transition that is not possible in sharp-interface and phase field models.

In this section we present solid-liquid phase change studies using model A, similar to those presented in section 3.3.3 for liquidsolid phase change. The space-time least squares formulation for a time strip (corresponding to an increment of time) with time marching is used to compute the evolution. Figure 17 shows a schematic of the space-time strip corresponding to the first increment of time, BCs and ICs, as well as dimensionless space-time domain and the dimensionless quantities.

0.9 0.8 0.7 Model A Constant ρ Variable ρ IC t=0.12 t=0.16 t=0.2

Model Problem 3: 1D Solid-Liquid Phase Change; Initiation and Propagation of Phase Transition

Figure 15: Model Problem 2: Evolution of thermal conductivity for liquid-solid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04 5

Figure 17: Solid-liquid phase transition: space-time strip, boundary conditions, and initial condition

Minimally conforming spaces are the same as described in section 3.3.3. Due to smoothness of the evolution, we choose k1 = 2 Figure 16: Model Problem 2: Interface location as a function of time, and k2 = 2 i.e. The of class C11 (Ω̄ext ), therefore the integrals in the C11 (Ω̄ext ), p = 9, ∆t = 0.04 STLSP are Lebesgue in x but Riemann in t. The space-time strip (∆t = 0.04) is discretized using 100 nine node space-time C11 (Ω̄ext ) 82

finite elements. Numerical studies were considered for the first space-time strip with phase change to determine adequate p-level for this discretization by starting with p-level of 3 (both in space and time) and incrementing it by two. At p-level of nine, I is of the order of 10−6 or lower and |(gi )|max ≤ 10−6 were achieved for all time steps. This ensures converged Newton’s linear method with line search as well as accurate evolution in the entire space-time domain. The numerical solutions computed using these values of h, p and k for [Ts , Tl ] = [−0.001, 0.001] are shown in figures 18–22. It may appear that presenting details of the evolutions of various quantities here is redundant in view of liquid-solid phase change model problem already considered, but this is not the case. In this case transition is from solid to liquid, thus evolutions of transport properties are quite different and hence essential to examine the resulting evolution of the solution.

Figures 18(a), 19(a), 20(a), 21(a), 22(a) show plots of T , L f , ρ, c p , and k versus x during the initial stages of the evolution (0 ≤ t ≤ 0.8). Continuous addition of heat from the right boundary progressively raises the temperature at the boundary and in the neighborhood of the boundary which eventually results in the initiation of phase change. Variations in L f (T ), c p (T ), k(T ) and ρ(T ) follow changes in temperature during evolution. From figure 18(a) we note that both constant and variable densities yield almost the same temperature distribution in the initial stages of the evolution. Figures 19(a), 21(a), and 22(a) show differences in the evolutions of L f , c p , and k for constant and variable density cases. As expected, variable density solutions lag constant density results, opposite of liquid-solid phase transition in section 3.3.3, model problem 2.

Figure 19: Model Problem 3: Evolution of latent heat for solid-liquid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figure 18: Model Problem 3: Evolution of temperature for solid-liquid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figures 18(b), 19(b), 20(b), 21(b), 22(b) show fully formed phase change transition region (solid to liquid) beginning with t = 3.2 and its propagation during evolution (3.2 ≤ t ≤ 19.2). For each space-time strip during time marching using ∆t = 0.04 ; I < O(10−6 ) and |(gi )|max ≤ 10−6 ensure accurate evolution that satisfies GDE quite well over the entire space-time domain of each space-time strip. All evolutions are smooth and free of oscillations. The influence of variable density can be seen more clearly in these graphs. The variable density evolution lags the constant density evolution for all values of time, and the difference between them increases as evolution proceeds. Here also we clearly observe linear heat conduction in the solid and liquid phases (constant but different slopes of T versus x) separated by a smooth transition region.

Similar to the liquid-solid studies presented in section 3.3.3, it is possible to use the solutions shown in figures 18–22 to follow the location of the phase transition front during the evolution. Figure 23 shows the location of the center of the transition zone for constant and variable densities. In contrast to similar results for liquid-solid phase transition shown in figure 16, here we note that the center of the phase transition zone for variable density case lags the constant density case during the entire evolution and the difference between the two increases as evolution proceeds. The phase transition for this model problem also cannot be simulated using sharp-interface and phase field models as this model problem requires initiation of phase transition that is not possible in sharp-interface and phase field models.

Figure 20: Model Problem 3: Evolution of density for solid-liquid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figure 21: Model Problem 3: Evolution of specific heat for solid-liquid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

In this section we consider liquid-solid and solid-liquid phase change in R2 using the mathematical model (47) (Model B).

In these numerical studies, we choose Model B, a system of first order PDEs that permits use of C0 local approximation in space and time. We consider a two dimensional domain in R2 consisting of a one unit square. A schematic of the domain, boundary conditions, and initial conditions are shown in figure 24. A constant heat flux is applied to each boundary (heat removal), except for the first time step in which heat flux changes continuously from zero at t = 0 to the constant value at t = ∆t.

x̂ =0.25 f t q̂y (x̂,0, tˆ) = −k̂ ∂∂T̂ŷ = −0.142 fBtu ;t ≥ ∆t t 2s

Figure 22: Model Problem 3: Evolution of thermal conductivity for solidliquid phase change, C11 (Ω̄ext ), p = 9, ∆t = 0.04

Figure 24: 2D liquid-solid phase transition: space-time slab, boundary conditions, and initial condition

A graded spatial discretization of the [1 × 1] spatial domain shown in figure 25 is constructed. Table 1 provides discretization 2 details of regions A, B, C and D. All four boundaries contain uni0 form heat flux q = −0.1 (cooling) for t ≥ ∆t. Evolution is com0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 puted (56 time steps) using p-level of 3 in space and time with Phase transition location, x ∆t = 0.0025 for the first 8 time steps and ∆t = 0.01 for the remainFigure 23: Model Problem 3: Interface location as a function of time, ing time steps. For this discretization, the C00 local approximation C11 (Ω̄ext ), p = 9, ∆t = 0.04 with p=3 yield I of O(10−6 ) or lower, confirming good accuracy of the solution. |gi |max ≤ 10−6 is used for convergence check in 4

yet. At t = 0.2 the entire boundary is in the transition zone except very small portions near y = 0 and y = 1 that have solidified. At t = 0.5 a significant portion of the boundary is completely frozen. Graphs of latent heat in figures 29(a) and (b) confirm these observations discussed here using figures 28(a) and (b). Graphs of the evolutions of ρ, c p , and k confirm these observations made from figures 28 and 29 and hence are not included. Evolutions are smooth and show that the differences between those with variable density and constant density are not as significant as for studies in R1 for the values of time reported here. As evolution proceeds, we expect more deviations between the two. As in the case of liquid-solid phase transition in R1 , here also the evolution with variable density leads the constant density evolution (more visible in figure 28(b) and 29). These studies demonstrate the strength of the work in moving front in R2 without front tracking techniques. In these numerical studies we have used [Ts , Tl ] = [−0.004, 0.004]. This model problem also cannot be simulated using phase field and sharp interface models due to the same reason as in the case of model problems in R1 . Symmetry of the evolution is quite obvious from figures 26 and 27.

the Newton’s linear method. For most time increments Newton’s linear method with line search converges in 5-10 iterations. Evolution of temperature T and latent heat L f calculated using variable density are shown in figures 26 and 27 using carpet plots for different values of time. Similar plots were generated for ρ, c p , k but are not shown for the sake of brevity. The carpet plots show evolutions to be oscillation free. Evolution and propagation of phase transition is demonstrated more clearly by using x, y plots of temperature and latent heat at the centerline and at the boundary. Figure 28(a) shows evolution of temperature at x = 0.5 (centerline) as a function of y for t = 0.01, 0.2, 0.5. Evolution of temperature T as a function of y at x = 0.0 (boundary) us shown in figure 28(b) for the same values of time. The evolution of latent heat L f for the same locations and for the same values of time are shown in figures 29(a) and (b). From figure 28(a) we observe that at t = 0.01, the phase transition has not initiated along the centerline. At t = 0.2, the portions of the domain closer to the boundary are experiencing phase transition. At t = 0.5, a significant length along y near the boundaries is in the transition zone with some portion near freezing. At the boundary, the situation is quite different (figure 28(b)). At t = 0.01 the phase transition has not initiated

Figure 26: Model Problem 4: Evolution of temperature for liquid-solid phase change in R2 , C00 (Ω̄xet ), p = 3, ∆t = 0.0025 for 0 ≤ t ≤ 0.02 and ∆t = 0.01 for t ≥ 0.02

Figure 27: Model Problem 4: Evolution of latent heat for liquid-solid phase change in R2 , C00 (Ω̄xet ), p = 3, ∆t = 0.0025 for 0 ≤ t ≤ 0.02 and ∆t = 0.01 for t ≥ 0.02

Figure 28: Model Problem 4: Evolution of temperature for liquid-solid phase change in R2 , C00 (Ω̄ext ), p = 3, ∆t = 0.0025 for 0 ≤ t ≤

0.02. and ∆t = 0.01 for t ≥ 0.02

Figure 29: Model Problem 4: Evolution of latent heat for liquid-solid phase change in R2 , C00 (Ω̄xet ), p = 3, ∆t = 0.0025 for 0 ≤ t ≤ 0.02 and ∆t = 0.01 for t ≥ 0.02

Model Problem 5: 2D Solid-Liquid Phase Change For this discretization, the C00 local approximations with p=3 yield I of O(10−6 ) or lower, confirming good accuracy of the solution. |gi |max ≤ 10−6 is used for convergence check of the Newton’s linear method. For most time increments Newton’s linear method with line search converges in 5-10 iterations. In these studies we have used [Ts , Tl ] = [−0.004, 0.004].

Here we also consider a two dimensional domain in R2 consisting of a one unit square. A schematic of the domain, boundary conditions, initial conditions and reference quantities are shown in figure 30. A constant heat flux is applied to each boundary, except for the first time step in which the heat flux changes continuously from zero at t = 0 to the constant value at t = ∆t. The graded discretization for the [1 × 1] spatial domain is same as in section 3.4.1, shown in figure 25, with details of regions A, B, C and D in Table 1. All four boundaries maintain uniform heat flux q = 0.1 (heating). Evolution is computed (50 time steps) using p-level of 3 in space and time with ∆t = 0.01.

Carpet plots similar to model problem 4 were also generated for this model problem with behaviors similar to model problem 4 and hence are not included here. Two dimensional line x, y plots of temperature T and latent heat L f are presented to demonstrate the phase transition more clearly.

constant L f further confirm completely liquid state of the matter.

Graphs of the evolutions ρ, c p , and k show evolutions that are in agreement with the evolutions of T and L f shown in figures 31 and 32 and hence are omitted for the sake of brevity.

In this case also the evolutions are smooth and show that the differences between the evolutions with variable and constant density are not as significant as for studies in R1 for the values of time reported here. As evolution proceeds we expect more deviations between the two evolutions. As in case of solid-liquid phase transition in R1 , here also the variable density evolution lags the constant density evolution (more visible in figures 31 and 32). This model problem also cannot be simulated using sharpinterface or phase field approaches as it requires initiation of phase transition.

Figure 30: 2D solid-liquid phase transition: space-time slab, boundary conditions, and initial condition

(a) Evolution of temperature at the centerline Figure 31(a) shows evolution of temperature at x = 0.5 (centerline) as a function of y for t = 0.01, 0.2, and 0.5. Evolution of 1 temperature T as a function of y at x = 0.0 (boundary) is shown in 0.9 figure 31(b) for the same values of time. The evolutions of latent heat L f for the same locations and for the same values of time are 0.8 shown in figures 32(a) and (b). From the evolution of temperature 0.7 in figure 31(a) we note that at t = 0.01, the phase transition has 0.6 not been initiated at the centerline. For t = 0.2 the entire region 0.5 Model B 0 ≤ y ≤ 1 is in the transition zone [Ts , 0]. At t = 0.5 the entire Constant ρ Variable ρ 0.4 zone 0 ≤ y ≤ 1 is still in the transition zone, but some portions near t=0.01 t=0.2 the boundaries are in [0, Tl ]. At the boundary (x = 0, 0 ≤ y ≤ 1) the 0.3 t=0.5 evolution of the temperature is quite different than at the centerline. 0.2 From figure 31(b) we find that at t = 0.01, the phase transition has 0.1 not commenced yet except in a small portion near y = 0 and y = 1 0 (horizontal boundaries at y = 0 and y = 1). At t = 0.2 the entire -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 length 0 ≤ y ≤ 1 is in the transition zone [0, Tl ]. At t = 0.5 a sigTemperature, T nificant portion of 0 ≤ y ≤ 1 near y = 0 and y = 1 is completely (b) Evolution of temperature at the boundary liquid. Graphs of latent heat L f in figures 32(a) and (b) confirm these Figure 31: Model Problem 5: Evolution of temperature for solid-liquid phase change in R2 , C00 (Ω̄xet ), p = 3, ∆t = 0.01 observations. In figure 32(b) we note that at time t = 0.5 the straight line portions of the graph near y = 0 and y = 1 meaning

given where x̄1 is the direction of the flow. If we choose x̄1 = x̄ and (0) (0) x̄2 = ȳ, v̄1 = ū, d σ̄x1x2 = d σ̄xy , then the dimensionless form of the mathematical model presented in section 72 (combined model) for this model problem reduces to      (0)  ∂ σ̄ ∂ p̄ τ p xy 0 0 d   − =0 f¯l  2 2  ∂ ȳ ρ0 v0 ∂ x̄ ρ0 v0         c̄ p (T̄ ) L f 0 ∂ L̄ f ∂ T̄  1 ∂ q̄y   ρ̄ + 2 + =0  Ec ReBr ∂ ȳ v0 ∂ T̄ ∂t (76)    
µ v ∂ ū  (0) 0 0   f¯l d σ̄xy = µ̄   L0 τ0 ∂ ȳ       ∂ T̄   q̄y = −k̄(T̄ ) ∂ ȳ

If we choose τ0 = ρ0 v20 , characteristic kinetic energy, then (76) reduces to  (0)  ∂ σ̄ ∂ p̄ xy  d  − =0 f¯l    ∂ x̄ ∂ ȳ        c̄ p (T̄ ) L f 0 ∂ L̄ f ∂ T̄ 1 ∂ q̄y    ρ̄ + 2 + =0  Ec ReBr ∂ ȳ v0 ∂ T̄ ∂t (78)  
µ̄ ∂ ū  (0)  f¯l d σ̄xy =    Re ∂ ȳ      ∂ T̄    q̄y = −k̄(T̄ ) ∂ ȳ

Figure 32: Model Problem 5: Evolution of latent heat for solid-liquid phase change in R2 , C00 (Ω̄xet ), p = 3, ∆t = 0.01

By substituting q̄y in the energy equation we can eliminate q̄y

3.5. Phase Transition Numerical Studies in the

as a dependent variable. We designate this as Model (a). Presence of Flow In this section we present numerical studies using mathematical model based on constant stress and zero velocity in the solid phase but nonzero velocity and stress field in the liquid and transition regions. The details of the mathematical model are presented in Section 2.3. In the following we present numerical results for fully developed flow between parallel plates in which the plates are being cooled to initiate and propagate liquid-solid phase transition. 3.5.1

Model Problem 6: Fully Developed Flow Between Parallel Plates

For this case we only need to consider evolution along any vertical line between the plates. The flow is pressure driven i.e. ∂∂x̄p̄ is 1

 (0)   ¯fl ∂ p̄ − ∂ d σ̄xy = 0    ∂ x̄ ∂ ȳ         c̄ p (T̄ ) L f 0 ∂ L̄ f ∂ T̄   ρ̄ + 2  Ec v0 ∂ T̄ ∂t (80)   2    1 ∂ k̄(T̄ ) ∂ T̄ ∂ 2 T̄   − + k̄(T̄ ) 2 = 0   ReBr ∂ ȳ ∂ ȳ ∂ T̄      
µ̄ ∂ ū  (0)  ¯fl d σ̄xy = Re ∂ ȳ

The mathematical model consists of (80) and (79) with ū, d σ̄xy , and T̄ as dependent variables. ∆tˆ = 12.5s , Model (b)

An alternate form of (80) can be derived by first substituting µ̄ ∂ ū (0) in the momentum equation and then recasting the d σ̄xy = Re ∂ ȳ momentum equation as a system of first order equations that enforce ∂∂ ūȳ = 0 in the solid region. We obtain the following:  (0) ∂ d τ̄xy ∂ p̄   ¯  Re fl − µ̄ =0   ∂ x̄ ∂ ȳ         c̄ p (T̄ ) L f 0 ∂ L̄ f ∂ T̄   ρ̄ + 2  Ec v0 ∂ T̄ ∂t (81)   2    ∂ 2 T̄ 1 ∂ k̄(T̄ ) ∂ T̄  − + k̄(T̄ ) 2 = 0    ReBr ∂ ȳ ∂ ȳ ∂ T̄      
∂ ū  (0)  ¯fl d τ̄xy = ∂ ȳ

Figures 33(a)-(c) show a schematic of the problem and spacetime strips for an increment of time ∆t from the lower plate to the center of the flow. Boundary conditions and initial condition are also shown in figures 33(b),(c). The lower plate is subjected to a temperature gradient of 0 to 0.3 (continuous and differentiable; This mathematical model consists of (81) and (79). In this cubic) for the first increment of time and held fixed thereafter as (0) shown in figure 33(d). Figure 33(e) shows a 40 element uniform ∂ τ̄ (0) model d τ̄xy , hence ∂∂ ūȳ and f¯l ∂∂ x̄p̄ = 0, and therefore d∂ ȳxy = 0 holds discretization for the space-time strip. Evolution is computed using in the solid region. This model is obviously an alternate way to the space-time discretization of figure 33(e) with time marching usachieve the desired physics of constant stress and zero velocity in ing solutions of class C11 in space and time with uniform p-levels the solid phase as in (80) and (79). of 9 in space and time. We consider both models (a) and (b) in the numerical calculaThe temperature range for the transition zone is chosen to be tions of the evolution. [T̄s , T̄l ] = [−0.003, 0.003]. ρ̄, c̄ p , k̄, and L̄ f are assumed to be con(0) ∂ τ̄ From these mathematical models it is clear that d∂ ȳxy = 0 in tinuous and differentiable functions of temperature in the transition (0) the solid region. This of course implies that d τ̄xy = C1 , a constant, zone. Fixed ∆t of 50.0 is considered during time marching. For com(0) in the solid region. We note that a constant d τ̄xy with f¯l = 0 and parison purposes, evolution is also computed using constant den∂ ū ∂ ȳ = 0 satisfies the last equation in the models (a) and (b) (i.e. last sity (ρ̄ = ρ̄l ). equation in (80) and (81)). Value of C1 is dictated by the deviatoric shear stress at the solid-liquid interface. Only when the deviatoric (0) shear stress at the liquid-solid interface becomes zero is d τ̄xy in Numerical Results Using Model (b) (0)

the solid region zero. In conclusion, a constant d τ̄xy in the solidified region is supported by the mathematical model. Its magnitude is largest at the initiation of freezing and is progressively reduced upon continued growth of the solidified region, eventually becoming zero when the entire flow domain is solidified. Transport Properties Once again, the solid phase is considered to be ice and the liquid phase is water, with the same properties that are listed in section 3.3.2. The reference and the dimensionless quantities are given as: ρ0 = ρ̂s

Figure 34 shows evolution of ū for t = 0 (IC), t = 1000, 2500, and 4000 for constant as well as variable density. Progressive increase in the solid zone that initiates at the plate is clearly observed as the evolution proceeds. With progressively increasing solid zone the flow height is progressively reduced. Since the flow is pressure driven ( ∂∂ x̄p̄ = −6.7182 × 10−5 , constant), this results in the progressive reduction in the flow rate. In other words, as evolution proceeds, the effective H2 is progressively reduced. (0)

Similar plots of temperature T̄ and d τ̄xy are shown in figures 35 and 36. In the solidified region as expected we observe linear heat conduction and constant deviatoric Cauchy shear stress of the same magnitude as at the solid-liquid interface. Figures 37-40 show evolutions of L̄ f , ρ̄, c̄ p , k̄ for variable as well as constant density for the same values of time. Smoothinterface approach in the transition region and the mathematical model work exceptionally well. 92

Figure 34: Model Problem 6: Evolution of velocity ū versus ȳ using model (b), C11 (Ω̄ext ), p = 9, ∆t = 50

0.3 Model (b) u– = 0, dσ–(0) xy = 0 in solid Constant ρ– Variable ρ– IC t = 1000 t = 2500 t = 4000

Figure 35: Model Problem 6: Evolution of temperature T̄ versus ȳ using model (b), C11 (Ω̄ext ), p = 9, ∆t = 50

∂ T̄ ∂ ȳ ȳ=0,t ∂ T̄ ∂ ȳ =0.3 0.5 Model (b) u– = 0, dσ–(0) xy = 0 in solid