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.

1. 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 Ç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 time-steps.

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.

2. Problem Formulations

2.1. Strong Formulation

The problem studied in this paper is based on the following strong formulation.

Find , , and defined on verifying the initial-boundary value problem defined by

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 .

2.2. 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 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 .

Assuming that , we consider the function such that   that is  , where is the trace operator.

For the existence of a solution, setting , we then seek satisfying the problem (1a)–(1e) with condition (1e) which can be considered now as homogeneous, according to the classical trace theory (see, e.g., [26]).

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.

Find  such that: where and where (according to (3), with notation for more convenience):

3. 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 to .
This is usually true in practice as shown in Figure 1 for an isothermal transformation. Under this assumption, there exists and such that , 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.

Proposition 1. The problem (9) is a particular case of the following model problem where (i)  and are two functional spaces; (ii)  is a bilinear form on ; (iii)  (where   denotes the dual space of ).

Proof. Let . Introducing the linear forms (we note the set of continuous linear map from to ):(i) such as ;(ii) such as ;the problem (9) is equivalent to the following problem

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.

Theorem 2 (Nečas). Let and two Hilbert spaces, , and . Then, the problem is well-posed if and only if:
Under those conditions, one get the estimation-inequality:

Proof. (see [28]).

According to our particular case, the Nečas theorem can be expressed as follows.

Corollary 3. Let and two Hilbert spaces, , and two bilinear forms.
Let .
The problem (9) is well-posed in the sense of Hadamard if and only if the two conditions below are satisfied:

Proof. Introducing the Banach operators associated to forms and , the problem (9) is equivalent to: and problem is well-posed if the conditions of the following theorem are verified [28].
Theorem 4  (see [28]).   is well-posed if and only if: (1)   is an isomorphism; (2)   is surjective,where   is the kernel of    and    is an operator such that   for all   and .
According to the lemma (see, e.g., [29] or [30]).
Lemma 5.  A Banach operator   is bijective if and only if: (i) , , ,   and(ii) , .
And like is reflexive, the two conditons (i) and (ii), that is the condition are equivalent to say that is an isomorphism.
Furthermore, considering the following lemma [29, 30].
Lemma 6.  Let . Then    is  surjective if and only if    such that  , .
And taking into account that is reflexive, 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  and of the Nečas theorem for the bilinear form can be more simply formulated as follows.

Theorem 7. 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 condition:

Proof. is a reflexive Banach space, and is coercive on (on ).
Thus, the self-adjoint Banach operator associated to is bijective, and then condition (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 , the Banach operator is clearly at least surjective, and satisfies in Theorem 7.

3.2. Existence and Uniqueness of for Fixed

Considering fixed, the following can be established.

Proposition 9. exists and is unique.

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: where the only independent variable is time .

3.3. 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 and according to the initial condition .

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 is a continuous mapping from to .

4. 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 to   .