4. Nonlinearity treatments

In order to solve the nonlinear problem we are considering numerically, it is necessary to deal with all the nonlinearities.

In our case, let us mention the strong non-linearity linked to the enthalpy function \(u(T)\) which takes into account the solid-liquid phase change, as well as the non-linearity linked to the possible presence of a condition at the limits of non-linear normal flow (radiation).

Recall that in the classical case of non-linear transient thermal problems without convection, i.e. \(V=0\), several resolution methods are proposed in the literature. There are methods using enthalpic formulations as well as methods using formulations at temperature, all of which aim to best treat the nonlinearity associated with enthalpy (phase change).

We refer the reader to reference [bib5] for a summary of the main methods encountered in the literature. However, it should be noted that due to the difficulty associated with the presence of the transport term \(V\mathrm{.}\mathit{grad}u(T)\) in the problem, none of these methods will be used later.

As in any iterative process, the aim of the numerical diagram in view is to find a temperature field \({T}^{n+1}\) at iteration \(n+1\), based on the temperature field \({T}^{n}\), the solution of the previous iteration.

4.1. Treatment of enthalpy-related nonlinearity

In order to deal with this nonlinearity, the strategy employed in this study was inspired by a technique for solving free border problems [bib3], which was originally proposed in [bib4].

Let’s consider enthalpy function \(u(T)\) as being given in a reciprocal form: Temperature as a function of enthalpy (inverse of function \(u(T)\)). In other words, we will have to deal with the following temperature-enthalpy relationship:

\(T=\tau (u)\) eq 4.1-1

The reason for this choice will be more clear in what follows. Indeed we will have to deal with a problem with two fields: a temperature field and an enthalpy field. The discretization of the inverse function [éq 4.1-1] makes it possible to increment the enthalpy field according to the current temperature field (and not the other way around) as follows:

The first-order expansion of function \(\tau (u)\) is as follows,

\({T}^{n+1}\mathrm{=}\tau ({u}^{n})+\tau \text{'}({u}^{n})({u}^{n+1}\mathrm{-}{u}^{n}),\) eq 4.1-2

where \(\tau \text{'}\) is the derivative of the function defined by [éq 4.1-1] in relation to its argument.

In order to take this nonlinearity into account, and starting with [éq 4.1-2], \({u}^{n+1}\) is replaced by an approximation as a function of the unknown temperature field \({T}^{n+1}\) in the following way:

\({u}^{n+1}\mathrm{-}{u}^{n}\mathrm{=}\omega ({T}^{n+1}\mathrm{-}\tau ({u}^{n})),\) eq 4.1-3

where \(\omega\) is a relaxation parameter, constant throughout the domain and throughout the iterative process, representing the term \(\frac{1}{\tau \text{'}({u}^{n})}\).

Because of the non-convexity of function \(\tau (u)\), this relaxation parameter must necessarily satisfy the following condition [bib2], [bib3]:

\(\omega \mathrm{\le }\frac{1}{\underset{n}{\text{max}}\tau \text{'}({u}^{n})}\) eq 4.1-4

In practice we take \(\omega \mathrm{=}\frac{1}{\underset{n}{\text{max}}\tau \text{'}({u}^{n})}\).

Taking into account the approximation [éq 4.1-3], the discretization of the second term of equation [éq3-2] is expressed in the following way:

\(\begin{array}{cc}\underset{\Omega }{\mathrm{\int }}V\text{.}\text{grad}{u}^{n+1}\psi d\Omega & \text{=}\underset{\Omega }{\mathrm{\int }}V\text{.}\text{grad}{u}^{n}\psi d\Omega +\underset{\Omega }{\mathrm{\int }}\omega V\text{.}\text{grad}{T}^{n+1}\psi d\Omega \mathrm{-}\underset{\Omega }{\mathrm{\int }}\omega V\text{.}\text{grad}\tau ({u}^{n})\psi d\Omega ,\end{array}\) eq 4.1-5

4.2. Treatment of nonlinearities related to the nonlinear Fourier condition and thermal conductivity

The nonlinearity related to the non-linear normal flow condition is treated by considering the first-order expansion of the function (assumed to be sufficiently regular) \(\alpha (T)\) which is given by:

\(\alpha ({T}^{n+1})\mathrm{=}\alpha ({T}^{n})+\alpha \text{'}({T}^{n})({T}^{n+1}\mathrm{-}{T}^{n}),\) eq 4.2-1

where \((\mathrm{.})\text{'}\) refers to the derivative of the \((\mathrm{.})\) function with respect to its argument.

It appeared necessary to decide on a strategy for discretizing the term \(k(T)\mathit{grad}T\) in equation [éq 3-2] in order to be able to deal with this nonlinearity for the stationary problem we are considering. To do this, we adopted the following approximation:

\(\begin{array}{cc}k({T}^{n+1})\text{grad}{T}^{n+1}& \text{=}k({T}^{n})\text{grad}{T}^{n+1}\mathrm{-}\left[k({T}^{n})\mathrm{-}k({T}^{n\mathrm{-}1})\right]\text{grad}{T}^{n}\end{array}\) eq 4.2-2

This discretization is in fact a simplification of the first-order development of the term \(k(T)\text{grad}T\). It is found to be effective in particular because of the low nonlinearity of function \(k(T)\) in practice.

Note:

Note also that the following purely explicit approximation:

\(k({T}^{n+1})\text{grad}{T}^{n+1}\mathrm{\approx }k({T}^{n})\text{grad}{T}^{n+1},\)

also gives satisfactory results. This observation was verified on the basis of several numerical experiments.