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

Abstract

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 Ω×[0,𝑇] verifying the initial-boundary value problem defined by𝑝=𝜆Δ𝑢inΩ,𝑢=𝑔(ℎ)inΩ,𝑝=𝜕ℎ𝜕𝑡inΩ,𝜆∇𝑢⋅𝑛=𝑗^(𝑝)onΓ_𝑗,𝑢=𝑢^(𝑝)onΓ_𝑢, (1a)(1b)(1c)(1d)(1e)

with the initial condition ℎ(𝑡=0)=ℎ₀. (2)

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 [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 𝑢^(𝑝).

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: ∀𝑣̃∈𝑉,_Ω𝑝𝑣̃d𝑉+_Ω𝜆∇𝑢⋅∇𝑣̃d𝑉−_Γ𝑗𝑗^(𝑝)𝑣̃d𝑆=0,∀𝑞∈𝑊,_Ω𝑢𝑞d𝑉−_Ω𝑔(ℎ)𝑞d𝑉=0, (3) with 𝑝=𝜕ℎ/𝜕𝑡 and where the two Hilbert spaces 𝑉 and 𝑊 are defined as 𝑉=𝐻¹_Γ𝑢(Ω)=𝑢∈𝐻¹(Ω),𝑢=0onΓ_𝑢,𝑊=𝐿²(Ω). (4)

Considering the above results, (3) can be rewritten in the following way.

Find (𝑢,𝑝,ℎ)∈𝑉×𝑊² such that: ∀𝑣∈𝑉,𝑎(𝑢,𝑣)+𝑏𝑣,𝑝=𝑙(𝑣),∀𝑞∈𝑊,𝑏𝑢,𝑞=𝑚ℎ,𝑞, (5)(6) where 𝑝=𝜕ℎ/𝜕𝑡 and where (according to (3), with notation 𝑢=𝑢 for more convenience): 𝑎(𝑢,𝑣)=_Ω𝑏=𝜆∇𝑢⋅∇𝑣d𝑉,𝑣,𝑝_Ω𝑙(𝑣)𝑝𝑣d𝑉,=_Γ𝑗𝑗^(𝑝)𝑏=𝑣d𝑆,𝑢,𝑞_Ω𝑚=𝑢𝑞d𝑉,ℎ,𝑞_Ω𝑔(ℎ)𝑞d𝑉. (7)

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 ℎ_fixed and 𝜅_ℎ such that 𝜅_ℎ=𝑔(ℎ_fixed), and this implies that the map 𝑚 in (6) is a unique linear form 𝑚𝑞=_Ω𝜅_ℎ𝑞d𝑉 (8) for ℎ fixed.
We consider in the sequel that the problem (5) and (6) is rewritten using (8) and becomes as follows.
Find (𝑢,𝑝)∈𝑉×𝑊 such that: ()(𝑣),𝑞.∀𝑣∈𝑉,𝑎𝑢,𝑣+𝑏𝑣,𝑝=𝑙∀𝑞∈𝑊,𝑏𝑢,𝑞=𝑚 (9)
To demonstrate the existence and uniqueness of (𝑢,𝑝), we first enunciate the preliminary result as follows.

Figure 1 

Typical evolution of the function 𝑔(ℎ) with recalescence.

Proposition 1. The problem (9) is a particular case of the following model problem 𝑞,𝐹𝑖𝑛𝑑𝑢∈𝐸𝑠𝑢𝑐ℎ𝑡ℎ𝑎𝑡∀𝑞∈𝐹,𝑘𝑢,𝑞=𝑓 𝑃_mod 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 ∀,.Find𝑢,𝑝∈𝐸=𝐹=𝑉×𝑊suchthat𝑣,𝑞∈𝑉×𝑊,𝑐𝑢,𝑝𝑣,𝑞=𝑑𝑣,𝑞 (10)

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 𝑃_mod is well-posed if and only if: ∃𝛼>0,inf_𝑣∈𝐸sup_𝑞∈𝐹𝑘𝑣,𝑞‖𝑣‖_𝐸𝑞_𝐹⩾𝛼,∀𝑞∈𝐹,∀𝑣∈𝐸,𝑘𝑣,𝑞=0⇒𝑞=0. 𝑁₁𝑁₂
Under those conditions, one get the estimation-inequality: ∀𝑓∈𝐹^′,‖𝑢‖_𝐸⩽1𝛼𝑓_𝐹′. 𝐸_𝐼

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 Ker(𝐵)={𝑣∈𝑉;∀𝑞∈𝑊,𝑏(𝑣,𝑞)=0}.
The problem (9) is well-posed in the sense of Hadamard if and only if the two conditions below are satisfied: ∃𝛼>0,inf_{(𝐵)𝑢∈Ker}sup_{(𝐵)𝑣∈Ker}𝑎(𝑢,𝑣)‖𝑢‖_𝑉‖𝑣‖_𝑉(𝐵),((𝐵)())⇒(),⩾𝛼,∀𝑣∈Ker∀𝑢∈Ker,𝑎𝑢,𝑣=0𝑣=0∃𝛽>0,inf_𝑞∈𝑊sup_𝑣∈𝑉𝑏𝑣,𝑞‖𝑣‖_𝑉𝑞_𝑊⩾𝛽. (Cond1)(Cond2)

Proof. Introducing the Banach operators associated to forms 𝑎 and 𝑏, the problem (9) is equivalent to: find𝑢,𝑝∈𝑉×𝑊suchthat𝐴𝑢+𝐵^𝑡𝑝=𝑙,𝐵𝑢=𝑚, (𝑃) 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)∃𝛼>0, ∀𝑣∈𝑉, ‖𝐴𝑣‖_𝑉⩾𝛼‖𝑣‖_𝑉,   and(ii)∀𝑣^′∈𝑉^′, (𝐴^𝑡𝑣^′=0)⇒(𝑣^′=0).
And like Ker(𝐵) is reflexive, the two conditons (i) and (ii), that is the condition (Cond1) 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  ∃𝛽>0   such that  ∀𝑞^′∈𝑊^′, ‖𝐵^𝑡𝑞^′‖_𝑉′⩾𝛽‖𝑞^′‖_𝑊′.
And taking into account that 𝑊 is reflexive, (Cond2) 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 (Cond2) condition: ∃𝛽>0,inf_𝑞∈𝑊sup_𝑣∈𝑉𝑏𝑣,𝑞‖𝑣‖_𝑉𝑞_𝑊⩾𝛽. (11)

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 (Cond1) (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 (Cond2) 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: 𝑝=𝜕ℎ()𝜕𝑡withinitialconditionℎ𝑡=0=ℎ₀, (12) 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 ℎ₀=ℎ(𝑥,𝑡=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 ℝ 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  ℝ.