Theoretical Framework of a Variational Formulation for Nonlinear Heat Transfer with Phase Changes

This paper discusses mathematical results for a variational formulation dedicated to heat transfer with phase changes. Practical finite element experiences show that the studied formulation can lead to difficulties for the numerical resolution at each time step. The aim of the paper is to show that such numerical pathologies do not come from the basic variational formulation by showing existence and uniqueness of the solution.


Introduction
Finite element analysis of heat transfer involving phase changes in solids is a numerical problem which has been extensively studied in the past years.The most classical example of the type of phenomena considered is the one of diffusion of heat in the presence of one or several phase transformations (fusion, solidification, metallurgical transformations of steels, etc.).A large body of literature has been devoted to this type of problems.Among the various approaches which have been proposed, three groups can roughly be distinguished: apparent heat capacity methods, fictitious heat source methods, and enthalpic methods [1].
Apparent heat capacity methods consist of artificially enhancing the heat capacity (derivative of the enthalpy with respect to the temperature) over the range of temperatures corresponding to the change of phase.This has been done for instance by Comini et al. [2], Morgan et al. [3], Lewis and Roberts [4], and C ¸etinel et al. [5], among many others.The application of the technique is straightforward when the phase change is spread over some relatively wide range of temperatures.Isothermal phase changes raise more problems since one may easily numerically miss the very high (virtually infinite) peak of heat capacity, thus failing to respect exact conservation of energy.One simple but costly remedy proposed by Rolph and Bathe [6] just consisted of using very small time-steps.Pham's [7] more elaborate proposal consisted of defining some "temperature correction" allowing to use relatively large time-steps while respecting strict conservation of energy.Pham's [8] comparison of various methods evidenced the high efficiency of this proposal.
Fictitious heat source methods consist of considering the heat absorbed or generated by the phase change as an internal heat source equal to the product of the latent heat of transformation and the transformation rate.This has been done by many authors, notably Hömberg [9], Murthy et al. [10], Han et al. [11], and Serajzadeh [12].Again, the application of this technique is simpler when the phase change occurs over a range of temperatures, but it can also deal with transformations occurring at constant temperature provided some precautions are taken.Voller [13] thus showed that the abrupt change of enthalpy resulting from the phase change can be accounted for accurately through explicit time-stepping or, with some appropriate linearization of the enthalpy over some arbitrarily narrow temperature range, implicit time-stepping.Another possibility proposed by Fachinotti et al. [14] consisted of accounting for the discontinuity of the enthalpy by identifying those elements where such a discontinuity occurs and performing exact, analytic integrations over them.
As explained by Feulvarch et al. [1], one trap with all this methods pertains to conservation of energy.In the most usual case, the temperature is a monotone function of time, so that the state of the phase change may be considered as a well-defined, single-valued function of temperature.In such a case, the apparent heat capacity and fictitious heat source methods are basically equivalent.Nevertheless, as noted by Pham [8], these methods are not equivalent for recalescence, which implies a non-monotone evolution of the temperature in time because the heating or the cooling is stopped.This phenomenon can only occur in the presence of some transformation kinetics; instantaneous transformations can only slow down but not stop the heating or the cooling.The equivalent heat capacity technique can be adapted to deal with such a case (see, e.g., Comini et al. [2] and Dalhuijsen and Segal [15]).However the fictitious heat source technique does the job in a more natural way since the only requirement is to know some appropriate expression of the transformation rate.The simulation of recalescence using this method started more than 20 years ago (see the review paper of Rappaz [16] and the works cited therein) and has become commonplace in the materials community since.It is generally based on explicit time-stepping, which prohibits the use of large timesteps.
Enthalpic methods consist of using the specific enthalpy instead of the temperature as principal nodal unknown.The heat flux is then expressed using the gradient of the enthalpy rather than that of the temperature, the temperature being considered as a function of the enthalpy.This technique dates back to Eyres et al. [17] and has been used since by Mundim and Fortes [18], Droux [19], Gremaud [20], and Pham [8], among others.It is particularly appropriate in the case of recalescence.Indeed, what makes the use of other methods cumbersome (though not impossible) in such a case is that the enthalpy becomes a multi-valued function of temperature, so that evaluating it from the temperature becomes a somewhat awkward task.On the other hand, use of the enthalpic technique raises no special difficulty since the temperature is always a well-defined, single-valued function of the enthalpy, so that deducing it from the enthalpy remains a straightforward task.One disadvantage of this technique, however, is that it implies a smooth spatial distribution of the enthalpy, which is unrealistic for transformations occurring at constant temperature.In materials with a sharp freezing point, for instance, a real, physical discontinuity exists on the freezing front.The enthalpic technique artificially spreads this discontinuity over one element length.A variant of the enthalpic method using both the enthalpy and the temperature as principal variables was proposed by Nedjar [21].
Although each method has shortcomings, many successful proposals have been made to numerically simulate heat diffusion with phase changes, even in the presence of recalescence with implicit time-stepping.To overcome this difficulty, a new method has been introduced by Feulvarch and Bergheau [22].The proposed formulation is based on the classical heat equation coupled with a function providing the temperature in terms of enthalpy.This enthalpy-temperature relation characterizes the kinetic of phase transformation and includes the latent heat.This approach has been extended to diffusion and precipitation of chemical elements by Feulvarch et al. [1] on the basis of mathematical results given by Leblond [23].This new implicit finite element technique rises sometimes some difficulties and a numerical solution can not always be found.The objective of this paper is to show that the basic variational formulation is well posed and leads to the existence and the uniqueness of a solution.Having established this result, one can conclude that the encountered difficulties are related to the numerical aspects.Find , , and ℎ defined on Ω×[0, ] verifying the initialboundary value problem defined by

Problem Formulations
with the initial condition In these equations,  is the temperature,  is the time derivative of the enthalpy ℎ;  is the thermal conductivity;  is the function providing the temperature  in terms of ℎ (equations (1b) cannot always be inverted to yield the enthalpy ℎ as a function of the temperature);  is the unit outward normal vector to the boundary;  () is a prescribed "input flux";  () is an imposed value of the temperature .

Weak Formulation.
In the problem (1a)-(1e), the variables ,  of space and time play different roles and this motivates a classical simulation approach which consists of a semi-discretization in space using finite elements which leads to solve ordinary differential equations in time, where the only independent variable is time  [24,25].Thus, we solve problem (1a)-(1e) by viewing (, ), (, ), and ℎ(, ) as functions defined on [0; +∞[ with values in a space , where  is a space of functions depending only on .
The weak formulation in space is obtained by multiplying (1a) by a weighting function , (1b) by a weighting function , and integrating over the domain Ω.
Since the boundary condition (1e) on Γ  is an inhomogeneous Dirichlet condition, let us introduce now an appropriate extension of  () .
Thereby, the weak formulation in space can be written as follows: with  = ℎ/ and where the two Hilbert spaces  and  are defined as Considering the above results, (3) can be rewritten in the following way.

Existence and Uniqueness of the Solution
To prove that the problem ( 5) and ( 6) leads to the existence and the uniqueness of , , and ℎ, we proceed in two steps, considering first a stationary spatial problem for  and , and then a time integration for ℎ.

3.1.
Existence and Uniqueness of (,) for ℎ Fixed.Considering ℎ is fixed, we place ourselves under the following working hypothesis.
Hypothesis 1.The function  in (1b) is a continuous mapping from R to R. This is usually true in practice as shown in Figure 1 for an isothermal transformation.Under this assumption, there exists ℎ fixed and  ℎ such that  ℎ = (ℎ fixed ), and this implies that the map  in ( 6) is a unique linear form for ℎ fixed.We consider in the sequel that the problem ( 5) and ( 6) is rewritten using (8) and becomes as follows.
Find (, ) ∈  ×  such that: To demonstrate the existence and uniqueness of (, ), we first enunciate the preliminary result as follows.
Thus, the problem ( 9) is well-posed in the sense of Hadamard if and only if the bilinear form  satisfies the two conditions of the Nečas [27] theorem.
According to our particular case, the Nečas theorem can be expressed as follows.and problem () is well-posed if the conditions of the following theorem are verified [28].
And like Ker() is reflexive, the two conditons (i) and (ii), that is the condition (Cond 1) are equivalent to say that  is an isomorphism.Furthermore, considering the following lemma [29,30].
Proof. is a reflexive Banach space, and  is coercive on  (on Ker()).Thus, the self-adjoint Banach operator  associated to  is bijective, and then condition (Cond 1) (which means  surjective and   injective) is verified.
Referring to the previous results, we can conclude with the following.Theorem 8.The previously stated problem (9) is well-posed, and admits a unique solution.
Proof.Since the bilinear form (, ) is coercive, and like the near form (, ) defines a norm in  =  2 (Ω), the Banach operator  is clearly at least surjective, and  satisfies (Cond 2) in Theorem 7. Proof.According to Theorem 8,  exists and is unique.From (1c) and ( 2), it can be easily stated that ℎ exists and is the unique solution of the temporal ordinary differential equation:

Existence and
where the only independent variable is time .

Existence and
Uniqueness of the Solution (,,ℎ) in Space and Time.Under Hypothesis 1, we showed in Section 3.1 that for ℎ fixed, we obtain a unique solution  which is the time derivative of ℎ at any point  of the domain Ω, independently of time  (cf.Theorem 8).
Thus, at any point  ∈ Ω, ℎ is the unique solution of the Cauchy problem (12) of Section 3.2 at each time t and according to the initial condition ℎ 0 = ℎ (,  = 0).
All these considerations allow us to state the following result to conclude this part.
Theorem 10.The previously stated problem (5) and (6) admits a unique solution (, , ℎ) in space and time, once function  in form (ℎ, ) = ∫ Ω (ℎ)d is a continuous mapping from R to R.

Conclusion
Mathematical results have been discussed for a variational formulation which is based on the heat equation coupled with a function providing the temperature in terms of enthalpy.Existence and uniqueness of the solution are established for this formulation which allows to model physical problems which cannot be easily computed with implicit time integration techniques.Therefore, we can conclude that computational difficulties encountered at each time step for the finite element method do not come from the basic variational formulation.In addition, all results are independent of the shape of the function  which may be nonlinear or not, since it is sufficient that  is a continuous mapping from R to R.

2. 1 .
Strong Formulation.The problem studied in this paper is based on the following strong formulation.
6. Let  ∈ L(, ).Then  is surjective if and onlyif ∃ > 0 such that ∀  ∈   , ‖    ‖   ⩾ ‖  ‖   .And taking into account that  is reflexive, (Cond 2) is equivalent to state that  is surjective.Moreover, one can easily demonstrate that the bilinear form  on  is coercive (it is a direct consequence of the Poincaré inequality).In such a case, conditons ( 1 ) and ( 2 ) of the Nečas theorem for the bilinear form  can be more simply formulated as follows.We postulate that the bilinear form  is coercive.Then, the problem (9) is well-posed if and only if the bilinear form  satisfies the previous (Cond 2) condition: