<figure><img src="image_1.png" /></figure>

HAL
open science

# The heat equation strictly written in enthalpy

Pierre Despret, Jean-Luc Dulong, Pierre Villon

► To cite this version:

Pierre Despret, Jean-Luc Dulong, Pierre Villon. The heat equation strictly written in enthalpy. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01514295

HAL Id: hal-01514295

https://hal.science/hal-01514295

Submitted on 26 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# The heat equation strictly written in enthlapy.

P. Despret<fnref n="1" />,<fnref n="2" />, J-L Dulong<fnref n="1" />, P. Villon<fnref n="1" />

**Résumé** — We propose here to modify the numerical enthalpy method developed by Feulvarch in 2009 in order to apply it to the solidification of alloys, with non linear material properties as a function of the temperature. We studied in the present paper the thermal conduction with phase change. We solve the energy balance equation written strictly in enthalpy. Usual numerical schemes are used. We display results compared with an analytical solution to prove the well convergence of the Modified Feulvarch's Method (MFM).

**Mots clés** — enthalpy method, phase change, non linear material properties, alloys

## 1. Introduction

To solve the phase transformation problem, notably in foundry process, E.Feulvarch [1] proposed to write a system of two equations with two unknowns the enthalpy and the temperature (h, T) as shown by the equation 1 to model the heat equation. He solved this system by an usual Newton Raphson Method (NRM).

$$ \begin{cases} \rho \frac{\partial h}{\partial t} = \nabla \cdot (k \nabla T) \\ T = g(h) \end{cases} \qquad (1) $$

The material is characterized by its thermal conductivity k, its density ρ and the function T=g(h) which is the bijection of the enthalpy curve as a function of the temperature (see 1). We can define the g function as following:

$$ h(x,t) = \hbar(T(x,t)) \Leftrightarrow T(x,t) = \hbar^{-1}(h(x,t)) = g(h(x,t)) \qquad (2) $$

ħ is the enthalpy depending of the temperature, given by the h-T curve. That means that the specific heat cₚ and the latent heat L are hidden inside the g function [2].

The Feulvarch's formulation has the main advantage, from our point of view, that the first equation becomes linear, with constant material properties. The enthalpy gap due to the release of the latent heat (L) is taken into account by the second equation.

## 2. Method

For our application with alloys and non linear material properties, the first equation becomes necessarily non-linear. So, we are going to use the previous formulation where the second equation is directly integrated into the first one, to obtain the so called Modified Feulvarch's Method (MFM) :

$$ \rho \frac{\partial h(x,t)}{\partial t} = \nabla \cdot [k(h(x,t)) \nabla g(h(x,t))] \qquad (3) $$

The non linearity is rejected on the right hand side, in the conduction term. Indeed, to restrict our phase transformation problem to a heat transfer equation and, so, to avoid the motion due to

<footnote n="1">¹Sorbonne universités, Université de technologie de Compiègne, UMR CNRS 7337, laboratoire Roberval, CR Royallieu CS60319 60203 Compiègne cedex, {pierre.despret, jean-luc.dulong, pierre.villon}@utc.fr, *corresponding author</footnote>

<footnote n="2">²MONTUPET SA, Service R&D, rue de Nogent, 60290 Laigneville</footnote>
contraction through the density variation, the density is supposed to be constant. All other material properties depend on the temperature, and so, through $g^{-1}(T)$, they can be calculated from the enthalpy $h$. We obtain, for the AlSi7G alloy, the following curves from the software JmatPro v8.

Moreover, this new formulation, depending from only one unknown, can be the starting point for a reduction of model resolution (PGD), developed in another paper.

The problem we are going to study in this paper is defined on $\Omega \times [0, t_{end}]$, where $\Omega$ is the space domain considered, $\Gamma$ the edge of this domain and $[0, t_{end}]$ the time of the simulation. It must verify the following boundary and initial conditions:

$$ 
\begin{aligned}
-k(h(\mathbf{x},t))\nabla g(h(\mathbf{x},t))\mathbf{n}=\varphi_N & \quad \text{if } \mathbf{x}\in\Gamma_N, \forall t \\
-k(h(\mathbf{x},t))\nabla g(h(\mathbf{x},t))\mathbf{n}=h_c(g(h(\mathbf{x},t))-T_\infty) & \quad \text{if } \mathbf{x}\in\Gamma_C, \forall t \\
h(\mathbf{x},t)=g^{-1}(T_D) & \quad \text{if } \mathbf{x}\in\Gamma_D, \forall t \\
h(\mathbf{x},0)=g^{-1}(T_0(\mathbf{x})) & \quad \text{at } t=0\text{s}
\end{aligned} 
\quad (4) $$

where $\varphi_N$ is the prescribed flux applied to the external boundary $\Gamma_N$, $h_c$ is the convective coefficient, $T_\infty$ is the temperature far from $\Omega$ applied in Robin condition on $\Gamma_C$. $T_D$ is the imposed Dirichlet temperature on $\Gamma_D$, with $\Gamma=\Gamma_N \cup \Gamma_C \cup \Gamma_C$ and $\mathbf{n}$ is the normal to the surface $\Gamma$.

The weak formulation gives the usual formulation:

$$ \int_{\Omega} \dot{h} h^* dV = - \int_{\Gamma_N} \varphi_N \cdot h^* dS - \int_{\Gamma_C} h_c(g(h) - T_\infty) \cdot h^* dS - \int_{\Omega} k(h) \cdot \nabla g(h) \cdot \nabla h^* dV \quad (5) $$

The spatial discretization is done by the standard Finite Element Method (FEM) where $N(\mathbf{x})$ is the shape function on each node:

$$ h(\mathbf{x},t) = \sum_{i} N_{i}(\mathbf{x}) h_{i}(t) = N^{T} \mathbf{h} \quad (6) $$

The function $g$ and its gradient, are discretized using the property of interpolation of finite element shape functions:

$$ 
\begin{aligned}
g(h(\mathbf{x},t)) &= g\left(\sum_{i} N_{i}(\mathbf{x}) h_{i}(t)\right) \approx \sum_{i} N_{i}(\mathbf{x}) g(h_{i}(t)) = \sum_{i} N_{i}(\mathbf{x}) g_{i}(t) = N^{T} g \\
\nabla g(h) &= \nabla\left(\sum_{i} N_{i} g_{i}\right) \approx \sum_{i} g_{i} \cdot \nabla N_{i}(\mathbf{x}) = [\nabla N] g
\end{aligned} 
\quad (7) $$

We define the usual matrices mass $[M]$ and conductivity $[K]$. The discretized heat equation is then, written strictly in enthalpy:

$$ [M]h + [K(h)]g(h) = F \quad (8) $$

Its time discretization is done by an implicit backward Euler scheme and it is linearized by the NRM procedure: for its $k^{th}$ iteration, the tangent matrix $A_T$ and the enthalpy increment $\delta h$ from the equation residual $r$ write:

$$ 
\begin{aligned}
[A_T(h_k)] &= -\frac{\partial r(h_k)}{\partial h} = [M] + \Delta t [K(h_k)] \frac{\partial g}{\partial h} + \Delta t \left[ \frac{\partial K}{\partial h} g(h_k) \right] \\
\delta h &= [A_T]^{-1} r(h_k)
\end{aligned} 
\quad (9) $$

At this point we must specify that we will not develop the non linearity from materials depending of temperature because it is done via the classic way. It must be noticed that for thermal/enthalpy material dependencies, for example the thermal conductivity, $k(H)$ has to be precalculated. We will
more develop the non linearity coming from g(H). Let's first define the [L] matrix whose i, j terms are:

$$[L]_{i,j} = \begin{bmatrix} \frac{\partial g}{\partial h} \\ \frac{\partial g_i}{\partial h_j} \end{bmatrix} = \frac{\partial g_i}{\partial h_j} = \frac{\partial g(h_i)}{\partial h_j} \quad \text{So that} \quad [L]_{i,j} = \begin{bmatrix} g'(h_i) & \text{if } i=j \\ 0 & \text{if } i \neq j \end{bmatrix} \qquad (10)$$

L matrix is diagonal whose each term is the slope of the g(H) function at each node i.

## 3. Results

The results of the previous described method will be displayed for an analytical solution. We do not compare MFM with the Feulvarch's one because both methods have the same root, so they lead to similar results. From these analytical equations, we used MATLAB® to solve the equations. To validate the MFM, we implemented it, and also the effective heat capacity method (solved with an implicit Newton-Raphson method), on a 1D Finite Element model.

<figure><img src="image_2.png" /><figcaption>FIGURE 1: SEMI INFINITE 1D BODY</figcaption></figure>

The model, described in Figure 1, is a semi infinite 1D solid [0; +∞] solved analytically in [3]. The material used is a simplified aluminium alloy. For material properties, please see the table 1 below. Due to the analytical model, the material properties are approximated by constants and the mass density value is the same whichever the temperature domain. Three domains of temperature are considered: the domain 1, below the solidus temperature ($T_2$ = 555°C), all is solid; the domain 2, defined in the mushy zone, that means between the solidus and the liquidus temperature ($T_1$ = 613°C) and the domain 3, above the liquidus temperature, all is liquid. With $L=446 000 J.kg^{-1}$. The Dirichlet temperature at x=0 is $T_D=0°C$, the initial temperature $T_0$ is uniform equal to 800°C. The solution is calculated on a 0.4m body length with one thermal measure point at x=0.08m.

<table><caption>Table 1: Summarize of the material properties</caption><tr><td>Liquid (Domain 1)</td><td>$k_1=180 W.m^{-1}.K^{-1}$</td><td>$\rho_1=2675 kg.m^3$</td><td>$c_{P 1}=900 J.K^{-1}$</td></tr><tr><td>Mushy(Domain 2)</td><td>$k_2=120 W.m^{-1}.K^{-1}$</td><td>$\rho_2=2675 kg.m^3$</td><td>$c_{P 2}=7433 J.K^{-1}$</td></tr><tr><td>Solid (Domain 3)</td><td>$k_3=85 W.m^{-1}.K^{-1}$</td><td>$\rho_3=2675 kg.m^3$</td><td>$c_{P 3}=1150 J.K^{-1}$</td></tr></table>

<figure><img src="image_3.png" /><figcaption>FIGURE 2: EVOLUTION OF THE TEMPERATURE IN FUNCTION OF TIME.</figcaption></figure>
This formulation implies a linear release of the latent heat between the liquidus and solidus temperatures. The two numerical methods (the usual effective heat capacity method $c_{peff}$ and the MFM) are compared with the Analytical method. We computed several simulations with the same material and boundary conditions but with different Maximum Time Step Allowed (MTSA). Indeed, the numerical simulations have been done with an adaptive time step. If the Newton-Raphson loop does not converge, the time step is divided by 4 and after 2 consecutive loops converged with less than 5 iterations, the time step is multiplied by 2, and so on, until to reach the MTSA.

The temperature evolution at $x=0.08m$, for a MTSA=1s, is plotted on the Figure 2. We can see the analytical method and the effective heat capacity method give very close results, superimposed on each other. The MFM follows the same path with a very small difference. To be more accurate in the comparison, we defined the global error $e$ as an $L_\infty$ norm in time and $L_2$ norm in space by equation 11. This allows us to evaluate the difference at the worst point.

$$e = \max_t \frac{\|T_{ana}(\mathbf{x},t) - T_{num}(\mathbf{x},t)\|_{L_2}}{\|T_{ana}(\mathbf{x},t)\|_{L_2}} \cdot 100 \quad (11)$$

In table 2, $e C_{peff}$ is the error as defined in the equation 11 between the effective heat capacity method and the analytical solution, while $e$ MFM is the error between the MFM and the analytical one. We also investigated the two numerical methods in terms of number of iterations. The gain of the MFM compared to the effective heat capacity is calculated by :

$$Gain = \frac{\text{number of iterations of the } c_{peff} \text{ method}}{\text{number of iterations of MFM}} \quad (12)$$

<table><caption>Table 2: Comparison of the gain and error for different MTSA.</caption><thead><tr><th>MTSA (s)</th><th>Gain</th><th>e C<sub>peff</sub> (%)</th><th>e Our method (%)</th></tr></thead><tbody><tr><td>0,50</td><td>3,12</td><td>0,08</td><td>0,98</td></tr><tr><td>1,00</td><td>10,82</td><td>0,09</td><td>0,97</td></tr><tr><td>5,00</td><td>47,46</td><td>0,25</td><td>0,94</td></tr><tr><td>10,00</td><td>30,13</td><td>0,25</td><td>0,93</td></tr><tr><td>50,00</td><td>3,37</td><td>0,42</td><td>1,16</td></tr></tbody></table>

This leads us to the table 2 where we can see the gain and the error as a function of the MTSA. If we look at the MTSA=1s, referring to the Figure 2, we can see that the error of the effective specific heat is very small regarding the analytical solution (0.09%) and the MFM is around 1%. If we look at the gain, the MFM is ten times faster. When we look globally the table 3, we are first attracted by the very high gains of the MFM, up to 47! It is true that the cost of the high speed computation is the error, around 1%. Please note that 1% acceptable as the computational methods are often above.

## 4. Conclusion

The MFM seems particularly promising to simulate phase transformation in material. This analytical solution allows us to justify that the MFM is first accurate, as seen on Figure 2, second robust and third efficient regarding the wildly used effective heat capacity method (table 2). The application of the MFM to an industrial tool is the topic of a coming publication.

## 5. References

[1] E. Feulvarch, An implicit finite elements algorithm for the simulation of diffusion with phase changes in solids, Int. J. Numer. Meth. Engng 2009; 78:1492-1512

[2] LEWIS and MORGAN, The finite element method in heat transfert analysis, John Wiley&sons, ISBN: 0471914240, chap V p 129

[3] CARSLAW and JEAGER, second Edition, Conduction of heat in solids, Oxford Science Publications, ISBN: 0198533683, Chap XI: Change of Phase p 282-286 et p290-291.