# Supercooling of phase change: A new modeling formulation using apparent specific heat capacity

Tanguy Davin, Benoit Lefez, Alain Guillet

► To cite this version:

Tanguy Davin, Benoit Lefez, Alain Guillet. Supercooling of phase change: A new modeling formulation using apparent specific heat capacity. International Journal of Thermal Sciences, 2020, 147, pp.106121. 10.1016/j.ijthermalsci.2019.106121. hal-02329539

HAL Id: hal-02329539

https://normandie-univ.hal.science/hal-02329539

Submitted on 20 Jul 2022

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.

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

Distributed under a Creative Commons Attribution - NonCommercial 4.0 International License
# Supercooling of phase change: A new modeling formulation using apparent specific heat capacity

Tanguy Davin, Benoit Lefez, Alain Guillet

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

## To cite this version:

Tanguy Davin, Benoit Lefez, Alain Guillet. Supercooling of phase change: A new modeling formulation using apparent specific heat capacity. International Journal of Thermal Sciences, 2020, 147, pp.106121. 10.1016/j.ijthermalsci.2019.106121. hal-02329539v2

HAL Id: hal-02329539

https://normandie-univ.hal.science/hal-02329539v2

Submitted on 16 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
1 Successful treatment of meningococcal bacteremia using oral doxycycline: a case report

2

3 S. Berbescu<sup>a,d</sup>, M. Matt<sup>a</sup>, A. Dinha<sup>a</sup>, S. Simeon<sup>a</sup>, L. Noussair<sup>b</sup>, L. Armand-Lefevre<sup>c</sup>, B. Davido<sup>a</sup>

4 <sup>a</sup>Infectious Diseases Unit, Raymond Poincaré Teaching Hospital, AP-HP, Versailles Saint-Quentin

5 University, Garches, France

6 <sup>b</sup>Microbiology Laboratory, Raymond Poincaré Teaching Hospital, AP-HP, Versailles Saint-

7 Quentin University, Garches, France

8 <sup>c</sup>Microbiology Laboratory, Bichat-Claude Bernard Teaching Hospital, AP-HP, Paris Diderot

9 University, Paris, France

10 <sup>d</sup>National Institute for Infectious Diseases “Prof. Dr. Matei Bals”, Carol Davila University,

11 Bucharest, Romania

12

13 <u>Corresponding author:</u>

14 B. Davido, Raymond Poincaré Teaching Hospital, Infectious Diseases, 104 Bd Raymond Poincaré,

15 Garches, 92380, France (benjamin.davido@aphp.fr), 01 47 10 77 74

16

17 <u>Keywords: meningococcal bacteremia, purpura fulminans, doxycycline</u>

18

19

20

21

22

23

24

25

26
# Supercooling of Phase Change: a new modeling formulation using apparent specific heat capacity

Tanguy Davin*, Benoît Lefez, Alain Guillet

tanguy.davin@univ-rouen.fr*

Groupe de Physique des matériaux, UMR CNRS 6634, Université de Rouen Normandie

5

## Abstract

The thermal behavior of Phase Change Materials (PCMs) is a major issue for cooling, heat storage and thermal management of various systems in general. Unfortunately many PCMs present supercooling which is a major drawback regarding the efficiency of cooling systems. Various solutions were proposed to model the thermal kinetics, in particular using the apparent specific heat $C_p(T)$ technique. But among them only few consider the supercooling effect. The present work considers this issue by focusing on the representation of the different steps of the supercooling phenomenon. This leads to different formulations of the apparent capacity $C_p(T, f_{super})$. The presented algorithm uses the lumped system analysis approach that is widely spread for first-order and multi-scale resolutions. The apparent capacity laws are explicitly presented for the different melting-crystallization steps. In particular, the punctually negative $C_p$ formulation, a helpful mathematical artifact, permits to reproduce the thermal dynamics during the local crystallization. The formulations and the models are discussed. Eventually, the temperature evolutions of an experimental system are compared to calculated data and the crystallization rate is compared to the literature.

10

15

20

**Keywords:** numerical simulation; lumped systems (capacitances) models; phase change; supercooling (subcooling); equivalent (apparent) specific heat capacity

## Nomenclature

<table><tr><td>a</td><td>Polynomial coefficients</td><td>-, 1/K or 1/K<sup>2</sup></td></tr><tr><td>A</td><td>Integrated specific heat area (enthalpy)</td><td>J/kg</td></tr><tr><td>[B]</td><td>Global heat sources vector</td><td>W</td></tr><tr><td>C</td><td>Heat capacitance</td><td>J/K</td></tr><tr><td>Cp</td><td>Specific heat capacity</td><td>J/(kg.K)</td></tr><tr><td>[D]</td><td>Global conductances matrix</td><td>W/K</td></tr><tr><td>E</td><td>Energy</td><td>J</td></tr><tr><td>erf</td><td>Error function</td><td>-</td></tr><tr><td>erf<sup>-1</sup></td><td>Inverse error function (erfinv)</td><td>-</td></tr><tr><td>f</td><td>Factor</td><td>-</td></tr><tr><td>G</td><td>Thermal conductance</td><td>W/K</td></tr><tr><td>h</td><td>Specific enthalpy</td><td>J/kg</td></tr><tr><td>H</td><td>Latent heat</td><td>J/kg</td></tr><tr><td>k</td><td>Thermal conductivity</td><td>W/(m.K)</td></tr></table>
27 Abstract

28 We report the case of an 18-year old immunocompetent man who presented to the hospital
29 with fever, headaches, and arthromyalgia that progressed to an erythematous rash, with a
30 history of a tick bite 72 hours earlier. The diagnosis of rickettsial infection was suspected and a
31 course of doxycycline was initiated for a total of 5 days. Evolution was rapidly favorable under
32 treatment, with resolution of symptoms within 24 hours. Blood cultures came back positive for
33 *Neisseria meningitidis* serotype B evocative of an authentic purpura fulminans.

34 Purpura fulminans is a medical emergency, a syndrome of intravascular thrombosis
35 characterized by a very rapid evolution that requires early recognition and specific treatment. It
36 is commonly described in the young and healthy patient and has high mortality and morbidity.
37 Common bacteria mainly associated with purpura fulminans are *Meningococcus spp.*,
38 *Pneumococcus spp.*, and *Staphylococcus spp.*

39

40 1. Case presentation

41 An 18-year old Caucasian male, belonging to the group of men having sex with men (MSM), was
42 admitted to the hospital due to fever (40°C), headaches and arthromyalgia that evolved with an
43 erythematous rash.

44 The recent patient history revealed that 2 weeks ago he visited the South of France (French
45 Riviera) with his parents, doing trips in the forest and collecting trash. Then he returned to Paris
46 to complete his high school graduation, followed by a couple of days of wild camping with
47 friends. While camping he was bitten by a tick on the back 72 hours before his admission.
<table><tr><td>N</td><td>Total number of nodes</td><td>-</td></tr><tr><td>Q</td><td>Heat sources</td><td>W</td></tr><tr><td>R</td><td>Thermal Resistance</td><td>K/W</td></tr><tr><td>t</td><td>Time</td><td>s</td></tr><tr><td>Δt</td><td>Time step</td><td>s</td></tr><tr><td>T</td><td>Temperature</td><td>°C</td></tr><tr><td>ΔT</td><td>Characteristic width of the Cp(T) Gaussian</td><td>°C</td></tr><tr><td>v</td><td>Speed rate</td><td>m/s</td></tr><tr><td>Z</td><td>1D-axis position</td><td>m</td></tr><tr><td>α</td><td>Convection coefficient</td><td>W/(m².K)</td></tr><tr><td>β</td><td>Factor of the negative Cp</td><td>-</td></tr><tr><td>ρ</td><td>Density</td><td>kg/m³</td></tr></table>

Abbreviations

<table><tr><td>ESHC</td><td>Enthalpy Specific Heat Capacity</td></tr><tr><td>PCM</td><td>Phase change Material</td></tr><tr><td>TC</td><td>Thermocouple</td></tr></table>

Subscripts

<table><tr><td>c</td><td>crystallization</td></tr><tr><td>eq</td><td>equivalent</td></tr><tr><td>ext</td><td>external</td></tr><tr><td>i,j</td><td>relative to the nodes</td></tr><tr><td>k</td><td>relative to the time step</td></tr><tr><td>l</td><td>liquid state</td></tr><tr><td>m</td><td>melting</td></tr><tr><td>mat</td><td>material</td></tr><tr><td>s</td><td>solid state</td></tr><tr><td>super</td><td>supercooling</td></tr><tr><td>tot</td><td>total</td></tr><tr><td>0</td><td>solid or liquid state</td></tr></table>

# 1 Introduction

Enthalpy of the liquid to solid (and vice-versa) phase change has been widely studied for thermal storage or cooling purpose [1]–[4]. PCMs (for Phase Change Materials) work as a thermal buffer, which load (on heating) and unload (on cooling) energy at an assumed constant temperature, in agreement with the Gibbs’ phase rule. Diverse numerical representations of phase change on various materials are already available in the literature [5]–[7]. In particular, two types of modeling have been provided by the literature, whether an Eulerian or Lagrangian approach is considered. One is based on a moving grid, permitting the interface tracking [8] and the other on fixed grid, introducing enthalpy methods. This former type of methods can be very accurate but it is relatively complicated especially for the discretization of the moving interface and it is time-consuming. These multi-domain methods involve the resolution of several phases. The latter considers a fixed discretization grid and only one set of equations for the different phases of the material. These one-domain methods, have been widely developed, as described in a recent exhaustive review on fixed-grid methods [9]. The idea of these methods is to represent the latent heat parameter in the formulation of the chosen variable. Indeed, the released or absorbed energy during the phase change is not located on a
The day after the tick bite he noticed an erythematous and purpuric rash (Figure 1) on his legs and pectoral area that was associated with fever and cephalalgia. Within 24 hours he presented an unusual fatigue, diarrhea, nausea and vomiting. At the Emergency Room, the clinical examination revealed a febrile patient, at 40°C, with hypotension (80/50 mmHg) and a generalized non pruritic macular rash that disappeared in vitropression. There were no neurological abnormalities.

Blood test showed elevated C reactive protein (150 mg/L) and neutrophilic leukocytosis, absence of renal dysfunction and normal transaminase values. Serologies for hepatitis C and B viruses, syphilis and HIV were negative.

In the hypothesis of a bacterial infection, antimicrobial therapy by 1 dose of amoxicillin-clavulanate of 875/125 mg was administered intravenously and the patient was transferred to our department without clinical improvement the day after.

The cerebrospinal fluid (CSF) was clear and colorless. The CSF cell count was 2/mm³ without red blood cells, hyperglycorachia (4.7 mmol/L) and normoproteinorrachia (0.22 G/L; reference value 0.20-0.40 G/L).

During the first hours of hospitalization, the appearance of the typical rash and the presence of tick bite history, made us suspicious of rickettsiosis and the patient got started on oral doxycycline 200 mg daily with rapidly a resolution of symptoms within 24 hours after his admission.

Blood culture at 48 hours came back positive for Neisseria meningitidis serotype B.
moving interface, but is included in the enthalpy variable of the system, either enthalpy-based (Enthalpy Methods EM), or Temperature-based (Enthalpy-Temperature Methods ETM) [10], [11]. The governing equations of enthalpy methods comprise the expression of the enthalpy variation by the mean of the specific capacity Cp(T) or enthalpy h(T).

Some authors directly use the enthalpy $h(T)$ [12], [13] while others prefer using the apparent specific capacity $C_p(T)$, leading to specific heat capacity (ESHC) models [14], also called zero latent heat models [14], [15]. Their implementation is hence simplified as the general equation remains unchanged. The actualization of the specific heat for every iteration is the only constraint. The apparent specific capacity variation permits to take into account the phase change transformation without modifying the general heat transfer equation and resolution. It is so convenient to consider this parameterization by using different formulation as Finite Element Method, Finite Difference or lumped systems, as it is proposed in the present study.

Supercooling also referred to as subcooling occurs when a compound is still liquid below the thermodynamic equilibrium temperature. This metastability implying a delayed phase change results from slow kinetics of crystallization. The nucleation/growth theory explains the formation of crystals[16]–[18]. It is possible to express the statistical degree of supercooling given some geometrical and microscopic properties as the homogeneous nucleation rate, the total number of atoms of the system, and the solid/liquid interfacial energy as well as the macroscopic latent heat and melting temperature [16], [19], [20]. For instance, for erythritol, a polyol which is used as an prototype of PCM with a severe supercooling in the experimental section of this article, the mean measured crystallization temperature for a 26g sample was 55°C: the degree of supercooling, difference between melting and crystallization temperatures is $T_m - T_c = 63°C$. The crystallization temperature $T_c$ was measured as 59°C with similar erythritol amount by Wei et al. [20], while the calculated homogeneous nucleation temperature, the lower limit of the crystallization temperature, is 46°C.

This theory also describes the crystal growth. The temperature gradients might be relatively high within a PCM block due to its poor conductivity (0.15-1 W/(m.K) usually [1], [2], [4], [21]). The growth of crystals can hence often be considered as heterogeneous. The moving interface can then be approximated in models, considering planar growth (1D solidification spread) as the present study although it is also possible to generalize to more complex configuration as proposed by Uzan [22]. Some laws indicating the crystallization rate as a function of the temperature and especially the degree of supercooling are interesting as it is needed for macroscopic representation.

Enthalpy methods can rely on an abundant literature for enthalpy-temperature characterization [23]. While the T-history method [24], [25] provides a simple identification of the enthalpy by comparing the free cooling of a sample with a reference, the Differential Scanning Calorimetry (DSC) is a precise tool to measure the enthalpy-temperature dependence. In particular, the control of the heat flux permits more flexibility and accuracy as both step and dynamic modes can be considered [26], [27]. Also, further studies on the DSC sample dynamics [28], [29], the identification of particular substances by inversion [30] or the comparison of different equipment have been achieved [31]. The accuracy of the method can be improved using a highly insulated device referred to as adiabatic scanning calorimeter [32]. Some authors [33], [34] have focused on supercooling using the T-history
An antimicrobial susceptibility testing was performed following the guidelines of the European Committee of Antimicrobial Susceptibility Testing (EUCAST) and showed that the strain was susceptible to all antibiotic tests, with the exception of a reduced susceptibility to amoxicillin that also applies to amoxicillin-clavulanate: MIC (minimum inhibitory concentration) > 0.5 µg/ml). In addition MIC to other antibiotics were as follows: cefotaxime (MIC <0.002 µg/ml), tetracycline (MIC 0.19 µg/ml), minocycline (minimum inhibitory concentration <0.064 µg/ml), tigecycline, ciprofloxacin and rifampicin.

Considering clinical and biological resolution under treatment (decrease from 18,000 WBC and a CRP level at 137 mg/L at day 0 to 4800 WBC and a CRP level at 10 mg/L at day 5), doxycycline was continued for a total of 5 days. Thereafter the patient was discharged with one dose of intravenous Ceftriaxone 1g administered for the eradication of pharyngeal carriage.

The patient was followed-up in consultation after one month for vaccination against other serotypes and further investigations revealed no complementary deficiency. He was considered cured with no sign of infection and normal biology. The patient declared the absence of other cases of infection in his close surroundings.

## 2. Discussion

In this case report, the atypical clinical presentation associated with the recent history of tick bite, was evocative of a potential rickettsial disease, considering the purpuric rash. So, it was decided to treat the patient with doxycycline. By luck, patient's strain was susceptible to cyclins which allowed full recovery.
methods. In particular, D'Avignon et al. [35] have compared the several ways commonly known of measuring the h(T) curves.

It is interesting to combine the crystal growth theory and the consideration of the enthalpy cooling curves during supercooling. This was done firstly by Günther et al. who proposed an enthalpy-based algorithm [34][12]. This study has two main purposes. The former is to investigate on the ways to quite easily implement the particular behavior of the phase change equation considering the formulation with apparent specific capacity h(T) (ESHC) instead of enthalpy h(T). The latter is the representativeness of the crystallization growth. Hence, the supercooling modeling is presented, discussed in terms of validity and robustness. Eventually, a comparison is done with experimental data.

# 2 Description of the model

## 2.1 Lumped system method

The lumped system analysis, also known as nodal method, is a relatively easy and adaptable formulation for the discretization and resolution. It is based on the resistance-capacitance electric analogy and allows to easily represent thermal links of a thermal system. The governing equations are described below. Based on the heat equation, the lumped system discretization leads to considered isothermal nodes. The main relationship taken into account is expressed as equation (1):

$$ \sum_{j} G_{ij}(T_{i} - T_{j}) + Q_{i} = C_{i} \frac{dT_{i}}{dt} \quad \forall i \in [1, N] \qquad (1) $$

Heat transfer between nodes is assured by the use of thermal conductances (inverse of resistances). These resistances represents the linear relation between heat flux and temperature difference of two nodes i and j. They can be used for conduction, convection or even radiation with a linearized relationship. The conductances are noted as Gij. The Qi term is gathering heat source and boundary conditions. The varying thermal mass i.e. heat capacitance Cj associated to every node guarantees the phase change. A typical lumped system network is represented in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Representation of a lumped system network for the discretization of the heat equation</figcaption></figure>

The equations system is processed with a matrix formulation implicit in space. For every time step all matrices are updated rendering the varying properties. Then the temperatures are calculated thanks to the time-explicit scheme in equation (2), where [D] is the global conductances matrix and [B] is the
Neisseria meningitidis is a gram-negative coccus acquired in the nasopharynx through
respiratory spread. The most common manifestations are related to meningitis or bacteremia,
but also it might result in localized infections such as arthritis or endocarditis (1,2).

One of the clinical manifestations of the meningococcal bacteremia is purpura fulminans, a rare
syndrome of intravascular thrombosis which occurs due to both infectious and non-infectious
causes. Infectious causes include meningococcal, streptococcal and staphylococcal sepsis, but
the most frequent is meningococcemia (3). Viral and non-infectious etiologies were also
described (primary varicella, HHV6 infection, deficiency of protein S or C).

The meningococcal disease is a life-threatening, but vaccine-preventable disease. The serogroups A, B, C, W, X and Y are responsible for the majority of meningococcal disease. In United States, the serogroup B caused 38% of all bacteremia included 69% of cases aged 18-23 years (4).

Using an alternative regimen to ceftriaxone, can be helpful in particular situations, notably in
patients allergic to beta-lactams. Even if ciprofloxacin is deemed to be effective, the resistance
of Neisseria group to fluoroquinolones was increased in the last decades, from 10% to almost
100% resistance, (5,6) particularly in men having sex with men (MSM). This led to other
therapeutic alternatives including aztreonam (7). However, aztreonam's spectrum activity does
not include gram-positive cocci such as *Streptococcus spp.* which is why using it can be very
dangerous when dealing with a purpura fulminans. Also, in case of unclear anamnesis,
meningococcal infection can be misdiagnosed such as illustrated in our case, and doxycycline
might find its place in such conditions.
110 global heat sources vector. The details on the mathematical transformation of equation (1) to (2) are given in the appendix.

$$ \left\{ \begin{array}{l} [\mathbf{D}] = [\mathbf{G}] - \frac{1}{\Delta t} [\mathbf{diag}(\mathbf{C}_i)] \\ [B] = \frac{1}{\Delta t} [\mathbf{diag}(\mathbf{C}_i)][T]^k + [Q] \\ [T]^{k+1} = -[\mathbf{D}]^{-1} \cdot [B] \end{array} \right. \qquad (2) $$

## 2.2 Specific capacity formulations for phase change without supercooling

The equations (2) are widely used for the calculation of lumped systems, but generally difficulties arise from the expression of the heat capacitances $C_i$. If the capacitances $C_i$ are simply the product of the node mass and specific heat capacity (following equation), the $Cp(T)$ expression is not simple

$$ C_i = m_i \cdot Cp(T_i) \quad \forall i \in [1, N] \qquad (3) $$

The latent energy of solidification can be added to the governing equations through the specific capacity $Cp$. For a pure substance as it is considered in the present study, the phase change occurs at a precisely reference melting temperature $T_m$. In theory, this would mean that the specific capacity $Cp(T)$ is a Heaviside function, represented in dotted black curve in Figure 2. The use of such a function would lead to huge numerical instabilities due to infinitely low time steps during phase change. A temperature range for the phase change can rule the problem, and is correct for mixture PCM, even if there is no thermodynamical justification to use such a mathematical formulation for pure substance. The simplest way to include the latent heat in the $Cp(T)$ equation is to add a constant within the $[T_m - \Delta T; T_m + \Delta T]$ range. This surface-equivalent square formulation is plotted in blue in Figure 2. To reduce errors and better represent the heat transfer (most enthalpy close to $T_m$), the Gaussian approximation is often preferred to the simple surface-equivalent square formula. The capacity is thus generally given by [12], [22], [36]:

$$ Cp(T) = Cp^{\text{gauss}}(T) = Cp^0 + \frac{H}{\Delta T \cdot \sqrt{\pi}} \cdot \exp\left(-\frac{(T - T_m)^2}{\Delta T^2}\right) \qquad (4) $$

The Gaussian standard deviation and latent heat are respectively expressed by $\Delta T$ and $H$. $\Delta T$ must be chosen as a compromise between the accuracy of the plateau phase change temperature and the numerical performance. Typically, a 1K deviation can be chosen for a reasonable numerical stability if this range is negligible compared to the temperature variation of the modeled system. For the present sketches, the range is deliberately enlarged. The integral on $T$ of the second term is equivalent to the latent heat $H$ as represented in Figure 2, also showing the simple square surface equivalence in blue. As the liquid and solid capacity can be different ($Cp(T)$ step at $T_m$), the regular specific capacity can be distinguished as $Cp^s$ or $Cp^l$. This is why there is a $Cp$ discontinuity (negligible compared to the maximum value) at $T_m$:

$$ \left\{ \begin{array}{l} Cp^0 = Cp^s \text{ in solid state (generally } T < T_m \\ Cp^0 = Cp^l \text{ in liquid state (generally } T > T_m \end{array} \right. \qquad (5) $$
110 Doxycycline has been proposed for decades in the prevention of bacterial meningococcal
111 meningitis and therefore is known to be effective against *Neisseria spp.* However, the molecule
112 is not usually tested in routine antibiograms, which in our case could be helpful.

113 Moreover, doxycycline in MSM is a common regimen well validated and prescribed against STIs,
114 especially in post-exposure prophylaxis (8). In the present case the patient was tested for
115 *Mycoplasma spp.* and *Neisseria spp.* even though the urinary sample was tested after he
116 received ceftriaxone and doxycycline, both came back negative. He was also negative for
117 syphilis and HIV.

118 Doxycycline is a lipophilic bacteriostatic antibiotic, which is completely absorbed after oral
119 administration but not recommended for the treatment of bacteremia, but could be considered
120 in an association. Our patient received an oral antibiotic, 2 times daily, for a total of 5 days.
121 Treating successfully bacteremia with bacteriostatic antibiotics was already described in
122 staphylococcal infections, but using Tigecycline, a derivate from doxycycline, in intravenous
123 administration (9).

124 Strains resistant to penicillin have been reported from a large number of countries and in
125 France the last report stated that 30% of *Neisseria meningitidis* are resistant to amoxicillin (10).
126 Overall our own experience describes an original alternative to treat Neisseria species.

127 The way to prevent morbidity and mortality following meningococcal disease is by vaccination.
128 Today the vaccinal recommendation includes the vaccine against A, C, W and Y serotype of
129 *Neisseria meningitidis*. Our patient was infected by a B serotype and currently there is still no
130 broad vaccine available covering all serotype B strains.
<figure><img src="image_4.png" /><figcaption>Figure 2: Common representation rendering the pure transformation, the Gaussian and the a simple surface-equivalent square approximations of the apparent specific capacity Cp(T) and the corresponding enthalpy h(T) without supercooling</figcaption></figure>

## 2.3 Specific capacity formulations for supercooling

The delicate formulation of the function Cp(T) with supercooling is the main purpose of this article and is described in the present section. In particular, two new sets of transition laws are proposed.

Günther et al. [12] proposed a method to represent the special variations of enthalpy during the supercooling process. In this study, the direct enthalpy $h(T)$ is used, as in a more recent study[37], [38]. It is interesting to study the different parameterization implied by the specific heat capacity and the laws that can be proposed to be implemented in numerous algorithm based on the direct Cp(T) formulation: all the ESHC models.

In the present work, the Gaussian approximation -relatively well spread in enthalpy based models- is considered as a basis to be modified during the cooling phase. In particular, the exit from the metastable state at the crystallization temperature $T_c$ leads to the definition of 2 different laws to reach back the stable state represented by the Gaussian approximation. These two different formulations of specific capacity lead to two equivalent phase change temperatures $T_{eq}^{(i \in \{1,2\})}$ in agreement with the conservation of the enthalpy. These equivalent temperatures represent the potential of the stored energy remaining after the heat transferred in the supercooled state (see Figure 4 and the following paragraphs). The different steps of the phase transitions with the supercooling are sketched in Figure 3. Among the two proposed sets of transition laws, the former is quite natural and corresponds to the simplest formulation: formulation 1 in blue, direct return to the stable state. The latter is more complex but might present some advantages: formulation 2 in green, for the representation of the crystallization dynamics.
131

132

### 3. Conclusions

133

134

135

136

137

This case report highlights the risk of misdiagnosis an authentic meningococcal disease presenting as a rickettsiosis. To the best of our knowledge this is the first case of meningococcal bacteremia treated successfully by oral doxycycline. Our case emphasizes that doxycycline could be an appropriate alternative therapy to intravenous ceftriaxone for particular non-severe situation or allergy to beta-lactams after susceptibility testing.

138

### Availability of data and materials:

139

140

All material and data described in the manuscript are available upon request to the corresponding author of the present article.

141

142

### Consent for publication:

143

144

Written informed consent was obtained from the patient for publication of this case report.

145

A copy of the written consent is available for review by the Editor of this journal.

146

147

### Ethics approval and consent to participate:

148

Not applicable.

149

150

Conflict of interest: none.

151
<figure><img src="image_5.png" /><figcaption>Figure 3: The different steps of the supercooling process implying the Cp(T) laws (formulation 1 in blue, formulation 2 in green)</figcaption></figure>

The principle of the supercooling modeling is to differentiate the specific capacity regarding the phase history. Due to the non-injectivity of the enthalpy function of the temperature, the complex phenomenon requires several laws with conditions. That is why a phase supercooling indicator $f_{super}$ must be defined. The common Gaussian law is used for the fusion always occurring at the same range of temperature around $T_m$ ($f_{super} = 0$). When liquid state is reached ($f_{super} = 1$), the specific capacity is assumed constant keeping the liquid property $Cp^l$. This metastable state is broken by the crystallization. This happens either if the crystallization temperature $T_c$ is reached or if the liquid-solid interface has arrived from the neighboring node. The tricky part is the following assumption on how to reach the equivalent temperature $T_{eq}$ (differentiated by a superscript in $T_{eq}^{(i)}$ presented in the following equations). The corresponding step $f_{super} = 2$ is represented by the arrows which can be more or less direct (blue or green paths respectively related to formulation 1 and 2). The proposed formulation 1 is found in the following equations (6). Also, the same Gaussian law is used for the end of cooling ($f_{super} = 0$). One can refer to the algorithmic conditions for the expression of $f_{super}$ in appendix.

$$ \left\{ \begin{array}{ll} \text{if } f_{super}(t) = 0 & Cp(T) = Cp^{gauss}(T) \\ \text{if } f_{super}(t) = 1 & Cp(T) = Cp^l \\ \text{if } f_{super}(t) = 2 & Cp(T) = Cp^{gauss}(T); T = T_{eq}^{(1)}; f_{super}(t + \Delta t) = 0 \end{array} \right. \qquad (6) $$

The formulation 1 is the simplest: the step $f_{super} = 2$ is skipped, and the temperature of the PCM node is directly set to the equivalent temperature $T_{eq}^{(1)}$. The stable Gaussian $Cp(T)$ curve is set for the next time iteration ($t + \Delta t$). In this case, the crystallization speed rate is assumed very high compared to the conduction rate. If this is relatively often true, some PCM like polyols present slow crystallization [39] that is interesting to represent. The only way to take the crystallization rate into account is to artificially limit the propagation of the crystallization. It can be done by a delay condition of a time $\Delta t_c = \Delta Z/v_c(T)$ on the supercooling factor for the neighboring node [12], [22], where $v_c$ is the crystallization rate and $\Delta Z$ is the distance between two PCM nodes.
152 Acknowledgments:

153 SB was in charge of conceptualization, data curation, methodology, project administration,
154 resources and writing – original draft; MM supervised the work; AD visualized the manuscript;
155 SS helped for the methods; LN gave resources for the case; LA helped for resources and
156 visualization; BD supervised, validated and wrote definite version of the manuscript.

157

158 4. References

159 1. GM Firestone. Meningococcus endocarditis. Am J Med Sci 1946; 211:556-564.

160 2. MA Apicella. Neisseria meningitidis GL Mandell, JE Bennett, R Dolin (Eds.), Principle
161 and practice of infectious diseases, Churchill Livingstone, Philadelphia (2000): 2228-
162 2241.

163 3. Kravitz GR, Dries DJ, Peterson ML et al. Purpura fulminans due to Staphylococcus
164 aureus. Clin Infect Dis 2005; 40:941-7.

165 4. Centers of disease Control and Prevention. Enhanced meningococcal disease
166 surveillance report, 2016: https://www.cdc.gov/meningococcal/downloads/NCIRD-
167 EMS-Report.pdf. (accessed on 2 October 2019).

168 5. Yang Y, Miao M, Gu W-M, et al. Antimicrobial susceptibility and molecular
determinants of quinolone resistance in *Neisseria gonorrhoeae* isolated from
169 Shanghai. Anti-microb Agents Chemother 2006; 58:868-72.

171 6. Goold PC, Bignell CJ. No way back for quinolones in the treatment of gonorrhea. Sex
172 Transm Infect 2006; 82:225-6.
<figure><img src="image_6.png" /><figcaption>Figure 4: Surfaces representing the equivalent enthalpy during heating and cooling (formulation 1)</figcaption></figure>

The equivalent phase change temperature $T_{eq}^{(1)}$ is calculated from the conservation of the enthalpy during the heating/cooling process. The total enthalpy can be represented by the surface integral of $Cp(T)$, plotted in Figure 4 for formulation 1. Here, $A_0$ corresponds to the enthalpy stored during heating from $T_1$ to $T_2$, two temperatures taken far from the melting temperature range. The total enthalpy restituted during cooling is the integral $A_1$ of the liquid or metastable state at $Cp^l$, plus the rest $A_2$ of the enthalpy in the form of the Gaussian law from $T_{eq}^{(1)}$. This condition gives the relationship to deduce the equivalent temperature $T_{eq}^{(1)}$. By subtracting surfaces, the simplified enthalpy area equivalence, represented in Figure 4 can be written as follows:

$$A_0 = A_1 + A_2 \Leftrightarrow A_0' = A_1'$$

$$Cp^l(T_{eq}^{(1)} - T_c) = \int_{T_{eq}^{(1)}}^{+\infty} \frac{H}{\Delta T \cdot \sqrt{\pi}} \cdot \exp\left(-\frac{(T - T_m)^2}{\Delta T^2}\right) \cdot dT \quad (7)$$

Using a change in variables and introducing the error function (calculation detailed in appendix), the following formula is obtained:

$$Cp^l(T_{eq}^{(1)} - T_c) = \frac{H}{2} \cdot \left( 1 - \operatorname{erf}\left(\frac{(T_{eq}^{(1)} - T_m)}{\Delta T}\right) \right) \quad (8)$$

Assuming $\Delta T$ is small enough so $Cp^l(T_{eq} - T_m)$ is negligible, the final expression of the equivalent temperature is, as also expressed by Uzan [22]:

$$T_{eq}^{(1)} = T_m + \Delta T \cdot \operatorname{erf}^{-1}\left(1 - \frac{2 \cdot Cp^l(T_m - T_c)}{H}\right) \quad (9)$$

The previous equation highlights the $\frac{Cp^l(T_m - T_c)}{H}$ term that corresponds to the ratio of the enthalpy transferred since the entry in the metastable state on the latent melting heat. This amount is representative of the solid fraction.

The idea of the second formulation is to intuitively reach the equivalent phase change temperature $T_{eq}^{(2)}$. This can be numerically done by setting a negative capacity during the step $f_{super}=2$ and so trying to represent the dynamics of crystallization. Though the idea of a negative $Cp$ may not be physically intuitive, some aspects of the interpretation of this negative $Cp$ are discussed in section 3.2.2. The formulation 2 is analogous to the first one but the law corresponding to $f_{super} = 2$ is different:
173 7. Davido B, Dinh A, Senard O, et al. Repurposing an old drug: aztreonam as a new treatment strategy for gonorrhea. J Antimicrob Chemother 2017; 72:1466-8.

175 8. Molina JM, Charreau I, Chidiac C, et al; Post-exposure prophylaxis with doxycycline to prevent sexually transmitted infections in men who have sex with men: Lancet Infect Dis 2018; 18:308-317.

178 9. Gardiner D, Dukart G, Cooper A, Babinchak T. Safety and efficacy of intravenous tigecycline in subjects with secondary bacteremia: pooled results from 8 phase. Clin Infect Dis 2010; 15:229-3.

181 10. Varon E. Epidemiology of acute bacterial meningitidis in adult patients in France. Med Mal Infect 2009; 39:432-44.

183

184

185
$$
\left\{
\begin{array}{ll}
\text{if } f_{\text{super}}(t) = 0 & Cp(T) = Cp^{\text{gauss}}(T) \\
\text{if } f_{\text{super}}(t) = 1 & Cp(T) = Cp^l \\
\text{if } f_{\text{super}}(t) = 2 & Cp(T) = -\beta \cdot Cp^l \quad ; \quad \text{if } T \ge T_{\text{eq}}^{(2)}, f_{\text{super}}(t + \Delta t) = 0
\end{array}
\right.
\qquad (10)
$$

Based on a new enthalpy surfaces represented in Figure 5, the enthalpy balance for formulation 2 can be done. In particular, the surface $A_2^*$ of the negative Cp is taken positive as the temperature is increasing: $A_0 = A_1 + A_2^* + A_2$. The relationship for the equivalent temperature in formulation 2 is given by equation (11), and $T_{eq}^{(2)}$ can be selectively expressed in equation (12) after the same assumption that $Cp^l(T_{eq} - T_m)$ is a negligible enthalpy:

$$
(1 + \beta) \cdot Cp^l (T_{eq}^{(2)} - T_c) = \frac{H}{2} \cdot \left( 1 - \operatorname{erf} \left( \frac{(T_{eq}^{(2)} - T_m)}{\Delta T} \right) \right) \qquad (11)
$$

$$
T_{eq}^{(2)} = T_m + \Delta T \cdot \operatorname{erf}^{-1} \left( 1 - \frac{2 \cdot Cp^l (1 + \beta) (T_m - T_c)}{H} \right) \quad (12)
$$

The speed of crystallization can be somehow modulated using a relative factor $\beta$. This factor can be
seen as the relative inverse of a crystallization rate. The lower value it takes, the quicker the
supercooled material reaches its phase change. For a very low factor $\beta$, the crystallization becomes
immediate and the formulation 2 is equivalent to the formulation 1. Noteworthy, for $\beta = 0$,
equations (9) and (12) are the same. The $\beta$ factor is discussed in more detail in section 3.2.2.

<figure><img src="image_7.png" /><figcaption>Figure 5: Surfaces representing the equivalent enthalpy (formulation 2)</figcaption></figure>

Noteworthy, for equation (12), a significant degree of supercooling and a high $\beta$ factor may lead to a numerical error value due to the use of the inverse function of the error function $erf^{-1}$ (erfinv in Matlab®) which has no input value outside [-1; 1]. Especially, when the sensible enthalpy (1 + $\beta$) $Cp^l(T_m - T_c)$ exceeds the latent enthalpy $H$, i.e. $T_{eq}^{(2)}$ is outside of the melting temperature rangea more simple relation can be preferred. Only on this case, equation (13) is used instead of (12):

$$
\text{if } (1 + \beta)Cp^l(T_m - T_c) > H \qquad T_{eq}^{(2)} = \frac{H + T_m \cdot (Cp^s - Cp^l) + T_c \cdot Cp^l \cdot (1 + \beta)}{Cp^l \cdot \beta + Cp^s} \quad (13)
$$

In practice, it may be better to use the full erf version of the equation to assess the equivalent
temperature. $T_{eq}^{(2)}$ is found using either the numerical iterative assessment of the solution of
equation (11), or selectively equations (12) and (13).
Enthalpy curves $h(T)$

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

Equivalent $C_p (T, f_{super})$ laws

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

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

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

Direct formulation

New formulation

1-node example

(both formulations)

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

Experimental/numerical comparison

(new formulation)

<figure><img src="image_7.png" /></figure>
Both formulations of the Cp(T) have been implemented and tested in the general algorithm (equation (2)). This implementation means equations (6) and (9) for formulation 1 and equations (10) and (12) for formulation 2. The two formulations have been compared for the case with only one node of PCM, then the relevancy of the results obtained with the formulation 2, as it is new, is deeply discussed in the following section.

## 3 Results and discussions

In this section, some results are given for a complete heating and cooling cycle. This permits to better understand the differentiation of the Cp(T) laws regarding the classic Gaussian law for a case without supercooling. As the lumped system models are often used as first order or multiscale approaches, the PCM one-node model is described in the separate first section 3.1. The particular supercooling laws can generate numerical problems especially for several nodes of phase change material. These problems are discussed in section 3.2.3.

The chosen PCM is erythritol as it is widely used because of its important energy density buffering (high *H*) but it also presents a high degree of supercooling and a relatively low crystallization rate. In order to compare the model with experimental measurements, a simple configuration is considered. An erythritol sample is heated from 50 to 150°C by a plate then cooled by the ambient external air at $T_{ext}=22°C$ through the same plate. The setup and its domain representation as thermal nodes are schematized in Figure 6. The heat transfer inside the PCM is considered one-dimensional so the solid-liquid interface is moving in a planar way. The PCM discretization node number in the orthogonal direction Z is an important parameter of the present study. Apart from the various PCM nodes due to its refinement, the rest of the lumped system network is unchanged. It is indeed composed of a few nodes representing the heating film, the diffusing plate, the containing capsule and the insulating parts. The elements in contact with the air are set with a Neumann boundary condition involving $T_{ext}$ and $\alpha_{ext}$, the natural convection coefficient. The materials as well as the insulation were taken into account through thermal resistances. Here only one node of PCM is drawn. For the PCM multi-nodes, the 6th node is replaced by a discretization of several nodes along Z-axis (slices of PCM in red dotted lines in Figure 6) as it is the temperature gradient direction. The thermophysical properties [39]–[41] and the boundary conditions are respectively listed in Table 1 and Table 2.

<figure><img src="image_8.png" /><figcaption>Figure 6: Sketch of the experimental test bench (left) and equivalent network for the nodal model (1 node of PCM)</figcaption></figure>

<table><thead><tr><td></td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Cp<sup>s</sup></td><td>1400</td><td>J/(kg.K)</td></tr><tr><td>Cp<sup>l</sup></td><td>2765</td><td>J/(kg.K)</td></tr><tr><td>H</td><td>333 000</td><td>J/kg</td></tr><tr><td>k<sub>s</sub></td><td>0.733</td><td>W/(m.K)</td></tr><tr><td>k<sub>l</sub></td><td>0.52</td><td>W/(m.K)</td></tr></tbody></table>
<table><caption>Table 1: Thermophysical properties values taken for Erytritol in the model</caption><tr><td>ρ</td><td>1480</td><td>kg/m³</td></tr><tr><td>T<sub>m</sub></td><td>118</td><td>°C</td></tr><tr><td>T<sub>c</sub></td><td>78</td><td>°C</td></tr></table>

<table><caption>Table 2: Boundary conditions and numerical parameters</caption><thead><tr><td></td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>ΔT</td><td>1</td><td>°C</td></tr><tr><td>Δt</td><td>0.25</td><td>s</td></tr><tr><td>T<sub>ext</sub></td><td>22</td><td>°C</td></tr><tr><td>α<sub>ext</sub></td><td>12</td><td>W/(m.K)</td></tr><tr><td>β</td><td>1 to 10</td><td>-</td></tr><tr><td>n<sub>PCM</sub></td><td>1 to 49</td><td>-</td></tr></tbody></table>

## 3.1 Case for a unique PCM node

Here are interesting results in the case of only one node of PCM. The whole block is then considered thermally thin and is represented by only one mean temperature. The model is stable as all tested parameters have led to convergent temperatures. Figure 7 represents the evolution of the variables of most interest. As they are expressed in different units and are of different orders of magnitude, one can refer to the specifications of the legend. Their evolutions are plotted as a function of the time during the cycle. The PCM mean temperature –at the only one node– and the plate temperature in contact with the PCM are given (node 3). It helps understanding the main heat transfer to the phase change material, by conduction. Also, the supercooling factor, the phase indicator and the apparent specific capacity are indicated.

<figure><img src="image_9.png" /><figcaption>Figure 7: Time graph representing the PCM and the contact plate temperatures, the apparent specific capacity and the phase indicators for both formulations (left: 1/ right: 2)</figcaption></figure>

In the heating period, the Gaussian apparent Cp leads to a relative plateau lasting for the whole melting process. This relative plateau (800 to 2000 s) is not completely flat since the latent heat is distributed over a temperature range and not exactly at $T_m$. The phase is considered as liquid as soon as the temperature reaches $T_m+2\Delta T$, (99.8% of the Gaussian area). The heating is continued a little
280

after the PCM is completely liquid. The liquid is then cooled till the crystallization temperature $T_c$ is reached ($t=3900$ s). Until then, both formulations are the same.

285

The direct formulation (#1) shows no transition to represent the crystallization dynamics. As the step $f_{super}=2$ is only one time step, the equivalent crystallization temperature $T_{eq}^{(1)}$ is instantly reached. At this crystallization time, more than half of the latent heat remains as can be seen on the black Cp curve and the relative plateau around $T_m$ is clearly visible (3900 to 4300 s). On the contrary the second formulation offers a gradual rise of the temperature from $T_c$ to reach the equivalent phase change temperature $T_{eq}^{(2)}$. Interestingly this transition when $f_{super}=2$ leads to a shortened relative plateau (4200 to 4400 s) as some of the energy is released during this period. As can be seen on the Cp curve, the part of the remaining latent heat is highly reduced (far less than the half) in this second case. The second formulation was tested using $\beta = 1$.

290

Both formulations must verify the conservation of energy. The balance involves the energy stored in all considered materials (internal), the energy given by the materials to the ambiance (boundary), and the energy source (production). The first two energies must be equal to the third ($E_{source} = E_{tot}$). These energies of the system are calculated for each time step. The energy stored in all materials $E_{mat}$ is defined as the mass-enthalpy product summed for all nodes, the PCM nodes taken separately considering the Gaussian law (from equation (4)) used most of the time:

295

$$E_{mat}(t) = \sum_{i \in PCM} C_i (T_i(t) - T_{t=0}) + \sum_{i \in PCM} m_i C_p g_{gauss}(T_i(t)). (T_i(t) - T_{t=0}) \quad (14)$$

300

The total energy is defined in equation (15) as the energy stored in the materials $E_{mat}$ plus the heat given to the ambiance (cooling power given to the boundary nodes). This term is calculated as the time-integral of the convection power. This energy is the sum for each node $i$ of the product of its convection conductance with the ambiance and the temperature difference with it. The energy balance is correct if the generated heat source $E_{source}$ is equal to this total energy.

$$E_{tot}(t) = E_{mat}(t) + \sum_{k} \sum_{i \_ convection} G_{i}^{amb} (T_{i}(t) - T_{amb}) \Delta t^{k} \quad (15)$$

305

In the left part of Figure 8 the energies $E_{mat}$, $E_{tot}$ and $E_{source}$ are plotted. The so-called source energy is the expected effective energy brought to the system, i.e. the integration of the step-like power condition. The break in both $E_{mat}$ and $E_{tot}$ curves corresponds to the ways in and out of the supercooled state.

310

The total energy balance is verified ($E_{source} = E_{tot}$) for both formulations except in the supercooled state. Indeed, in equation (14), the expression of the stored energy $E_{mat}$ considers the simple integration of the Cp using the law with the Gaussian peak at $T_m$ (equation (4)). Inside this metastable state, due to the complexity of the representation of the actual stored energy, it is not possible to verify the balance with this Gaussian law.

The right part of Figure 8 also presents the real Cp(T) evolutions during the whole cycle as it was sketched in Figure 3. The incoherence in the energy balance starts as soon as the explicit Cp diverges from the Gaussian law. It is important to notice that before the crystallization and as soon as it is finished, $E_{source} = E_{tot}$.
<figure><img src="image_10.png" /><figcaption>Figure 8: Energy balance verification (left) and explicit Cp (right) for both formulation</figcaption></figure>

## 3.2 PCM multi-nodes model

### 3.2.1 Behavior of the multi-nodes model

The PCM is then discretized to better represent the temperature gradient inside the material. The number of PCM nodes in the Z direction was studied. For stability purposes (see section 3.2.3) the $\beta$ factor could not be kept equal to 1 for such a large range of node numbers. In this example $\beta$ = 10. For an easy representable example, the variables are plotted in Figure 9 for a restricted number of nodes (7 elements). It helps understanding the numerical hypotheses and results, especially by interpreting the local Cp. First, it can be noted that the Gaussian profiles of the Cp(T) are distorted during the heating process. Even at the coldest side, the phase change is numerically initiated at an early stage of the relative plateau. Though, the interface, which can be defined as the location where the temperature is equal to $T_m$ (Cp peaks), moves quasi-linearly. More interestingly for our study, the recrystallization shows two types of supercooled node behaviors. At first, the local supercooling degree and the $\beta$ factor are so high that $T_{eq}^{(2)}$ is far from reaching the melting temperature range: $(1 + \beta)Cp^l(T_m - T_c)$ exceeds the latent enthalpy H (equation (13) is then used for the equivalent temperature). The heat released during the transition phase $f_{super}=2$ permits to increase in temperature the crystalizing node and the neighboring nodes by conduction. As the phase change occurs, the local crystallization temperature then increases and for the last three crystallizing nodes, the melting temperature range is reached (local temperature close to $T_m$), showing the rest of the Gaussian Cp.
<figure><img src="image_11.png" /><figcaption>Figure 9: Time graph representing the main variables for the 7-nodes case, formulation 2</figcaption></figure>

The influence of the nodes number was tested to see the influence on the temperature at the middle of the PCM. The energy balance is still verified, as shown by the right part of Figure 10. Here in the left part, the temperatures in the middle of the PCM and the plate are given. The inside temperature is more accurately calculated by increasing the number of nodes, until expectedly reaching a node-independent temperature evolution. In this case, a few dozens elements are enough for a reasonable precision; the 25 nodes mesh leads to insignificant error compared to the one with doubled elements (0.1% difference on $E_{tot}$ at the final state). It is also interesting to notice that the model is almost unchanged outside the phase change. Thereby, the domain must either focus on one node for simplicity or at least a reasonable amount of node for a relevant representation of phenomena during the phase change.

<figure><img src="image_12.png" /><figcaption>Figure 10: Influence of the number of nodes on the energy and the characteristic temperatures: (1-3-7-25-49) nodes, formulation 2</figcaption></figure>

The $C_p(T)$ laws must be considered more cautiously in order to get a good interpretation of the supercooling curves, especially for a refined mesh.
### 3.2.2 Physical interpretation of the apparent Cp

This section is the study of the meaning of the expression of the apparent specific heat capacity and its link to the macroscopic crystallization rate. Some experimental work have used the formulation of the apparent specific heat capacity to characterize the phase change [42]. Also using the Cp(T) formulation during supercooling, D'Avignon et al. [35] found that the explicit capacity experimentally identified could be negative during the cooling process. The data also show that the general shape of h(T) during the supercooling phenomenon seems similar to the following Figure 11 representing the present hypothesis (formulation 1 in blue and formulation 2 in green). In such studies representing h(T) [35] [43] [44] [45] the equivalent step 2 ($f_{super}=0$) of the supercooling process shows either null or slightly negative slope, respectively concordant with the formulation 1 and 2. Indeed, the derivation of these experimental enthalpy curves during supercooling present punctually negative Cp. In their study, D'Avignon et al. note the unrealistic nature of the negative Cp and the thermal transfers occurring inside the sample have to be considered.

Though it is not correct in classic thermodynamics, the idea of a negative Cp has been previously raised by astronomers observing that adding energy to a star or a star cluster would make it expand and cool down [46]. The negative specific heat has also been found in black holes [47]. The theory has then been developed to show that it could actually be verified at a microscopic level [48]–[50] in isolated or nearly-isolated systems. Though, these microscopic considerations could produce large fluctuations in non-equilibrium states as phase transition. In particular, it could lead to an abnormal slope of the supercooled or superheated matter corresponding to the negative specific heat [51]. Other researchers have noticed that negative Cp was commonly encountered in DSC during glass transition [52]–[55]. They have developed several-level models which consider this non-equilibrium state. In particular, negative Cp could occur when the system relaxation time to reach the internal energy barrier is long compared to the observation time.

The relevancy and the actual limits of these physical considerations are complex to determine and might be further studied. Especially, in the macroscopic scale as is presently dealt, the negative Cp should not be considered as more than a mathematical artifact. However, its value is not meaningless as it influences the crystallization dynamics. Indeed, in the present paper, the factor $\beta$ of the punctually negative Cp can be rewritten as the inverse of a relative crystallization rate factor $f_c$. $a_0$, $a_1$ and $a_2$ are the polynomial coefficients defining the crystallization factor dependency of temperature.

$$ \beta = \frac{1}{f_c} \quad (if \ f_{super} = 2) \qquad (16) $$

$$ f_c = a_0 + a_1(T_m - T_c) + a_2(T_m - T_c)^2 $$

Günther et al. [12] have used a second order polynomial of temperature formula for the assessment of the speed of solidification. This definition was also used here and was more precisely expressed as the local supercooling degree ($T_m - T_c^{local}$).
<figure><img src="image_13.png" /><figcaption>Figure 11: Representation of the apparent specific capacity Cp(T) with supercooling for formulation 1 (blue) and 2 (green)</figcaption></figure>

It is interesting to see the influence of this factor (inverse of the slope of the enthalpy curve, called the crystallization rate factor $f_c$) and verify the dependence with the actual modeled crystallization rate. Various laws for the crystallization factor $f_c$ have been tested and represented in Figure 12 for the cooling step. Here only the polynomial expression of $f_c$ has been changed, keeping all properties of erythritol. The polynomial coefficients are given in appendix. The resulting temperatures and the crystallization variables are plotted as functions of the local crystallization temperature: the implemented law $f_c(T_c^{local})$ and the actual crystallization rate $v_c$. The actual crystallization rate is the directly calculated as $\Delta Z / \Delta t_c$: the distance between two nodes divided by the time between the two nodes crystallizations (occurring at $T_c^{local}$). Thereby by setting different coefficients of the polynomial form of this crystallization rate factor, it is possible to retrieve several types of typical cooling curves given by Garg et al.[56] and Safari et al. [57]. Also, the results on the crystallization rate are interesting. It can be noted that the plotted theoretical crystallization rate is equal to $f_c$ affected by a constant coefficient in order to be represented homogeneously with the actual crystallization rate. The important observation is that both curves show the same trend, as expected for the correct representation of the cooling curve. The proportionality factor between the crystallization factor and the crystallization rate is affected by a lot of parameters, and firstly the cooling rate.

In practical, there might be two ways for parameterization. If the crystallization rate $v_c$ law is known, the crystallization factor might be assessed by using the same law, then adjusting the proportional factor by computing and checking the actual rate. If a representative enthalpy-temperature curve is known, the crystallization rate factor $f_c$ can be assessed by fitting the curve. The considered curve, obtained for example with DSC, must be taken in relevant conditions. That is to say in the same range of cooling rate as proposed by Solomon et al. [58].
<figure><img src="image_14.png" /><figcaption>Figure 12: Influence of the law of the supercooling factor $f_c(T)$ (formulation 2) on the temperature and the actual crystallization rate</figcaption></figure>

### 3.2.3 Limitations due to instabilities

The different Cp laws may lead to problems in robustness of the general algorithm regarding the chosen parameters and conditions. Especially, the transition with local and punctual high contributions of energy may produce instabilities. As warned by Korin et al. [59], the general algorithm of the supercooling model of Günther et al. [12] should comprise some verifications for stability purpose. First, the results of the equivalent temperatures with the error function or its inverse -erf or erfinv in Matlab®- (equations (8-9-11-12)) must be checked as well as the balance of the total energy (equation (15)). As described below for formulation 2, a negative capacity might generate instabilities. The main parameters found to influence the stability were the time step and the crystallization rate factor.

<figure><img src="image_15.png" /><figcaption>Figure 13: Instabilities (a) at the very beginning of the supercooling: improving with a smaller time step (b) of by forcing a lower $f_c$ (c), formulation 2</figcaption></figure>

As the resolution algorithm (equation (2)) uses a time-explicit scheme, it is also important to care about the time step. The beginning of the supercooling can be critical as shown in Figure 13. In some cases, the temperatures diverge. Actually the stability is given by the conditions which must be in agreement with the calculation characteristics. A smaller time step stabilizes the calculation for this particular moment. It should be chosen relevantly as it affects the calculation time. In fact, the
characteristic crystallization time of all supercooling nodes must be at least a few times bigger than the time step. That is why the algorithm can be improved in stability by decreasing the time step, but also parameterizing a lower crystallization rate. It might then be good to set a maximum value of the $f_c$ coefficient. It limits the very high rates but makes the algorithm more robust.

A very low $f_c$ coefficient can lead to errors in the energy balance. Also, in some case the crystallization can be obtained from two locations. In particular, for really slow crystallization, the crystallization temperature can be reached at a node different from the current solidifying node. For these reasons, it might be judicious to set a minimum crystallization factor $f_c$. As seen in Figure 12, a pretty wide range of crystallization rate can nevertheless be obtained without these restrictions.

## 4 Experimental validation

The experimental bench was composed of a heating film comporting a Joule effect circuit. This film was stuck to a diffusive aluminium plate in contact with the thin capsule containing the PCM sample. The aluminium capsule was 52mm wide in outside diameter and 0.4mm thick while the 5mm thick diffusing plate diameter was 90mm. The PCM sample is composed of 26g of erythritol provided by Nu3®, which is spread inside the whole capsule. With the density taken in the solid state, as it is assumed constant in the model, the PCM volume thus represents a cylinder with a 51.2mm diameter and a 12.6 mm height. The PCM capsule was insulated from the top so it was heated then naturally cooled through the plate. The heating plate temperature was regulated by an on/off regulator to generate cycles between 50 and 150 °C. The flat resistance delivered a power of 22W at 20°C, with a temperature resistivity dependence. The bench was instrumented with thermocouples: one stuck on the plate measuring the assumed position of node 3. Four other thermocouples were inserted inside the PCM at different heights in Z. Taking Z=0 at the capsule bottom, Z=1mm for TC2 and TC3 (at the center and decentered), Z=4.5mm for TC4 and Z=9mm for TC5 (close to the PCM free surface). Dozens of cycles were endured in order to see the repeatability of the crystallization. The cycles presented discrepancy due to the random nature of the crystallization (various $T_c$). The degree of supercooling was often more than 50 °C. A cycle with an average crystallization temperature (78°C) was chosen for the comparison.

The described elements (thermo-physical properties [60] and geometry) and the ambient air convection were taken into account in the capacitances model. In order to better represent the general heat transfer, an algorithm of inverse method was adapted to refine the air convection coefficient, some thermophysical properties and the coefficients of the polynomial expression of $f_c$. The direct model is the representation of the experimental setup with the formulation 2. The inverse algorithm is not detailed nor discussed in this article. The quadratic error made on one cycle is about 3 °C. This error seems reasonable since the assumptions (neglected internal convection, 1D in PCM, constant density) do not permit a very accurate representation. Indeed, the measurements during the heating step show that temperatures tend to reach together at the end of the melting phase, showing a mixing effect. This is corroborated by the overestimated value of the identified liquid conductivity (0.52 versus 0.326 W/(m.K) in the literature) and by the Rayleigh number, assessed around 5000 considering the biggest dimension of the PCM block and the literature [61].
<figure><img src="image_16.png" /><figcaption>Figure 14: Comparison of temperatures between experimental measurements and simulation (formulation 2 with identification “Id”) during cooling</figcaption></figure>

As it is shown in Figure 14, the identified modeled temperatures (thick lines) are in good agreement with the measurements (thin lines). In particular, it is important to identify the polynomial coefficients to well represent the heat transfer during the crystallization. The identification focused on the cooling shows a good behavior, maybe except the very beginning of the crystallization. At that point, the crystallization is very quick as the local degree of supercooling is high. After, for the upper nodes, the crystallization rate is highly reduced and the model better describes the thermal dynamics, in particular for the light blue and pink curves.

For measured configurations with various crystallization temperatures, the identification of the coefficients of the polynomial expression of $f_c(T)$ was made. Here the second order polynomial law (equation (16)) was used but some other crystallization factor laws might be tested for better fit. The resulting crystallization rates are plotted in Figure 15 with a reference law provided by Duquesne et al. [39] considering the growth of erythritol. This comparison highlights a reasonably similar behavior between the identification model and the literature.
<figure><img src="image_17.png" /><figcaption>Figure 15: Comparison of the crystallization rate between the present results with the identification model and the literature [39]</figcaption></figure>

# 5 Conclusions and prospects

This article presents a new formulation (referred to as formulation 2) for the apparent specific heat capacity model (ESHC) with supercooling. It is an enriching alternative to existing models (referred to as formulation 1) which do not directly represent the crystallization kinetics. Indeed, the punctually negative *Cp* is –though not explained by the classic thermodynamics– permits to consider the temperature growth in a recrystallizing material returning latent energy to the system. This mathematical artifact also enables to introduce the crystallization rate.

This new formulation 2 can be implemented in any apparent specific heat capacity (ESHC) model. In the present study, it has been associated to the lumped system approach. Both coarse and refined PCM domains have been tested. With only one node, the first order approach often used with lumped systems, can yet represent bulk effects and guarantees a good stability. The multi-nodes case reveals instability issues, but also some possibilities with the crystallization rate laws through the factor $f_c(T)$.

The multi-nodes model has been compared to an experimental assumed 1D system where the inner temperatures dynamics were measured. The low crystallization rate being particularly adequate for the proposed formulation, erythritol, a widely studied PCM, was chosen as for experimental comparison. The calculations are in good agreement with the measured temperatures and the crystallization rate found in the literature.

# Acknowledgments

The presented work has been conducted for the CRIOS project. This project was supported by the Normandy region and the European Union. Europe is committed in Normandy through the FEDER funding (Fond Européens de Développement Régional).
515

References

520

[1] F. Agyenim, N. Hewitt, P. Eames, and M. Smyth, "A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS)," *Renew. Sustain. Energy Rev.*, vol. 14, no. 2, pp. 615-628, Feb. 2010.

[2] N. R. Jankowski and F. P. McCluskey, "A review of phase change materials for vehicle component thermal buffering," *Appl. Energy*, vol. 113, pp. 1525-1561, Jan. 2014.

[3] S. K. Sahoo, M. K. Das, and P. Rath, "Application of TCE-PCM based heat sinks for cooling of electronic components: A review," *Renew. Sustain. Energy Rev.*, vol. 59, pp. 550-582, Jun. 2016.

525

[4] B. Zalba, J. M. Marín, L. F. Cabeza, and H. Mehling, "Review on thermal energy storage with phase change: materials, heat transfer analysis and applications," *Appl. Therm. Eng.*, vol. 23, no. 3, pp. 251-283, Feb. 2003.

[5] V. Alexiades, *Mathematical modeling of melting and freezing processes*. CRC Press, 1992.

[6] T. A. Kowalewski and D. Gobin, *Phase change with convection: modelling and validation*. Springer, 2004.

530

[7] N. S. Dhaidan and J. M. Khodadadi, "Melting and convection of phase change materials in different shape containers: A review," *Renew. Sustain. Energy Rev.*, vol. 43, pp. 449-477, 2015.

[8] J. M. Hyman, "Numerical methods for tracking interfaces," *Phys. Nonlinear Phenom.*, vol. 12, no. 1-3, pp. 396-407, 1984.

535

[9] A. König-Haagen, E. Franquet, E. Pernot, and D. Brüggemann, "A comprehensive benchmark of fixed-grid methods for the modeling of melting," *Int. J. Therm. Sci.*, vol. 118, pp. 69-103, 2017.

[10] V. R. Voller, C. R. Swaminathan, and B. G. Thomas, "Fixed grid techniques for phase change problems: a review," *Int. J. Numer. Methods Eng.*, vol. 30, no. 4, pp. 875-898, 1990.

[11] C. R. Swaminathan and V. R. Voller, "On the enthalpy method," *Int. J. Numer. Methods Heat Fluid Flow*, vol. 3, no. 3, pp. 233-244, 1993.

540

[12] E. Günther, H. Mehling, and S. Hiebler, "Modeling of subcooling and solidification of phase change materials," *Model. Simul. Mater. Sci. Eng.*, vol. 15, no. 8, p. 879, 2007.

[13] C. Le Bot and D. Delaunay, "Rapid solidification of indium: Modeling subcooling," *Mater. Charact.*, vol. 59, no. 5, pp. 519-527, May 2008.

545

[14] K. A. Fikiin, "Generalized numerical modelling of unsteady heat transfer during cooling and freezing using an improved enthalpy method and quasi-one-dimensional formulation," *Int. J. Refrig.*, vol. 19, no. 2, pp. 132-140, 1996.

[15] Y. Jarny and D. Delaunay, "Numerical resolution of a phase change problem with zero latent heat," *Numer. Heat Transf. Part B Fundam.*, vol. 16, no. 1, pp. 125-141, 1990.

550

[16] D. Turnbull and J. C. Fisher, "Rate of nucleation in condensed systems," *J. Chem. Phys.*, vol. 17, no. 1, pp. 71-73, 1949.

[17] I. V. Markov, *Crystal growth for beginners: fundamentals of nucleation, crystal growth, and epitaxy*. World Scientific, 2003.

555

[18] M. Qian, P. Cao, M. A. Easton, S. D. McDonald, and D. H. StJohn, "An analytical model for constitutional supercooling-driven grain formation and grain size prediction," *Acta Mater.*, vol. 58, no. 9, pp. 3262-3270, 2010.

[19] D. Turnbull, "Formation of crystal nuclei in liquid metals," *J. Appl. Phys.*, vol. 21, no. 10, pp. 1022-1028, 1950.

[20] L. Wei and K. Ohsasa, "Supercooling and solidification behavior of phase change material," *ISIJ Int.*, vol. 50, no. 9, pp. 1265-1269, 2010.

560

[21] A. Sharma, V. V. Tyagi, C. R. Chen, and D. Buddhi, "Review on thermal energy storage with phase change materials and applications," *Renew. Sustain. Energy Rev.*, vol. 13, no. 2, pp. 318-345, Feb. 2009.

[22] A. Y. Uzan, Y. Kozak, Y. Korin, I. Harary, H. Mehling, and G. Ziskind, "A novel multi-dimensional model for solidification process with supercooling," *Int. J. Heat Mass Transf.*, vol. 106, no. Supplement C, pp. 91-102, Mar. 2017.

565
[23] A. Ristić et al., "Engineering and processing of PCMs, TCMs and sorption materials," *Energy Procedia*, vol. 91, pp. 207-217, 2016.

[24] Z. Yinping and J. Yi, "A simple method, the-history method, of determining the heat of fusion, specific heat and thermal conductivity of phase-change materials," *Meas. Sci. Technol.*, vol. 10, no. 3, p. 201, 1999.

[25] A. Lázaro, E. Günther, H. Mehling, S. Hiebler, J. M. Marin, and B. Zalba, "Verification of a T-history installation to measure enthalpy versus temperature curves of phase change materials," *Meas. Sci. Technol.*, vol. 17, no. 8, p. 2168, 2006.

[26] C. Barreneche, A. Solé, L. Miró, I. Martorell, A. I. Fernández, and L. F. Cabeza, "Study on differential scanning calorimetry analysis with two operation modes and organic and inorganic phase change material (PCM)," *Thermochim. Acta*, vol. 553, pp. 23-26, 2013.

[27] C. Castellón, E. Günther, H. Mehling, S. Hiebler, and L. F. Cabeza, "Determination of the enthalpy of PCM as a function of temperature using a heat-flux DSC—A study of different measurement procedures and their accuracy," *Int. J. Energy Res.*, vol. 32, no. 13, pp. 1258-1265, 2008.

[28] A. Greco and A. Maffezzoli, "Correction of melting peaks of different PE grades accounting for heat transfer in DSC samples," *Polym. Test.*, vol. 27, no. 1, pp. 61-74, 2008.

[29] G. Albright, M. Farid, and S. Al-Hallaj, "Development of a model for compensating the influence of temperature gradients within the sample on DSC-results on phase change materials," *J. Therm. Anal. Calorim.*, vol. 101, no. 3, pp. 1155-1160, May 2010.

[30] E. Franquet, S. Gibout, J.-P. Bédécarrats, D. Haillot, and J.-P. Dumas, "Inverse method for the identification of the enthalpy of phase change materials from calorimetry experiments," *Thermochim. Acta*, vol. 546, pp. 61-80, 2012.

[31] A. Lazaro et al., "Intercomparative tests on phase change materials characterisation with differential scanning calorimeter," *Appl. Energy*, vol. 109, pp. 415-420, 2013.

[32] J. Leys, P. Losada-Pérez, E. Slenders, C. Glorieux, and J. Thoen, "Investigation of the melting behavior of the reference materials biphenyl and phenyl salicylate by a new type adiabatic scanning calorimeter," *Thermochim. Acta*, vol. 582, pp. 68-76, 2014.

[33] B. Sandnes, "The physics and the chemistry of the heat pad," *Am. J. Phys.*, vol. 76, no. 6, pp. 546-550, 2008.

[34] E. Günther, S. Hiebler, H. Mehling, and R. Redlich, "Enthalpy of Phase Change Materials as a Function of Temperature: Required Accuracy and Suitable Measurement Methods," *Int. J. Thermophys.*, vol. 30, no. 4, pp. 1257-1269, Aug. 2009.

[35] K. D'Avignon and M. Kummert, "Assessment of T-History Method Variants to Obtain Enthalpy-Temperature Curves for Phase Change Materials With Significant Subcooling," *J. Therm. Sci. Eng. Appl.*, vol. 7, no. 4, p. 041015, 2015.

[36] K. Darkwa and P. W. O'Callaghan, "Simulation of phase change drywalls in a passive solar building," *Appl. Therm. Eng.*, vol. 26, no. 8, pp. 853-858, Jun. 2006.

[37] A. Yehya, H. Naji, and L. Zalewski, "Simulation numérique de la convection avec surfusion dans un matériau à changement de phase via la méthode de Boltzmann," *SFT Conf.*, no. 122, 2016.

[38] A. Yehya and H. Naji, "Thermal Lattice Boltzmann Simulation of Entropy Generation within a Square Enclosure for Sensible and Latent Heat Transfers," *Appl. Sci.*, vol. 5, no. 4, pp. 1904-1921, 2015.

[39] M. Duquesne, A. Godin, E. P. del Barrio, and F. Achchaq, "Crystal growth kinetics of sugar alcohols as phase change materials for thermal energy storage," *Energy Procedia*, vol. 139, pp. 315-321, 2017.

[40] R. A. Talja and Y. H. Roos, "Phase and state transition effects on dielectric, mechanical, and thermal properties of polyols," *Thermochim. Acta*, vol. 380, no. 2, pp. 109-121, 2001.

[41] A. Kaizawa et al., "Thermophysical and heat transfer properties of phase change material candidate for waste heat transportation system," *Heat Mass Transf.*, vol. 44, no. 7, pp. 763-769, May 2008.
[42] E. D. Kravvaritis, K. A. Antonopoulos, and C. Tzivanidis, "Improvements to the measurement of the thermal properties of phase change materials," *Meas. Sci. Technol.*, vol. 21, no. 4, p. 045103, 2010.

620

[43] L. Desgrosseilliers, P. Allred, D. Groulx, and M. A. White, "Determination of enthalpy-temperature-composition relations in incongruent-melting phase change materials," *Appl. Therm. Eng.*, vol. 61, no. 2, pp. 193-197, Nov. 2013.

[44] B. Sandnes and J. Rekstad, "Supercooling salt hydrates: Stored enthalpy as a function of temperature," *Sol. Energy*, vol. 80, no. 5, pp. 616-625, May 2006.

625

[45] H. P. Garg, S. C. Mullick, and A. K. Bhargava, "Latent Heat or Phase Change Thermal Energy Storage," in *Solar Thermal Energy Storage*, Springer, Dordrecht, 1985, pp. 154-291.

[46] A. S. Eddington and A. S. Eddington, *The internal constitution of the stars*. Cambridge University Press, 1988.

630

[47] S. W. Hawking, "Black holes and thermodynamics," *Phys. Rev. D*, vol. 13, no. 2, pp. 191-197, Jan. 1976.

[48] D. J. Wales and J. P. K. Doye, "Coexistence and phase separation in clusters: From the small to the not-so-small regime," *J. Chem. Phys.*, vol. 103, no. 8, pp. 3061-3070, 1995.

[49] D. Lynden-Bell, "Negative specific heat in astronomy, physics and chemistry," *Phys. Stat. Mech. Its Appl.*, vol. 263, no. 1-4, pp. 293-304, 1999.

635

[50] R. M. Lynden-Bell and D. J. Wales, "Free energy barriers to melting in atomic clusters," *J. Chem. Phys.*, vol. 101, no. 2, pp. 1460-1476, 1994.

[51] D. Lynden-Bell and R. M. Lynden-Bell, "On the negative specific heat paradox," *Mon. Not. R. Astron. Soc.*, vol. 181, no. 3, pp. 405-419, 1977.

640

[52] J. Bisquert, "Master equation approach to the non-equilibrium negative specific heat at the glass transition," *Am. J. Phys.*, vol. 73, no. 8, pp. 735-741, 2005.

[53] M. A. DeBolt, A. J. Easteal, P. B. MACEDO, and C. T. Moynihan, "Analysis of structural relaxation in glass using rate heating data," *J. Am. Ceram. Soc.*, vol. 59, no. 1-2, pp. 16-21, 1976.

[54] R. K. P. Zia, E. L. Pr\a estgaard, and O. G. Mouritsen, "Getting more from pushing less: Negative specific heat and conductivity in nonequilibrium steady states," *Am. J. Phys.*, vol. 70, no. 4, pp. 384-392, 2002.

645

[55] T. Tao, T. Odagaki, and A. Yoshimori, "Cooling rate dependence of specific heat in systems out of equilibrium," *J. Chem. Phys.*, vol. 122, no. 4, p. 044505, 2005.

[56] H. P. Garg, S. C. Mullick, and A. K. Bhargava, "Latent Heat or Phase Change Thermal Energy Storage," in *Solar Thermal Energy Storage*, H. P. Garg, S. C. Mullick, and A. K. Bhargava, Eds. Dordrecht: Springer Netherlands, 1985, pp. 154-291 (section 3.2, p.160-170).

650

[57] A. Safari, R. Saidur, F. A. Sulaiman, Y. Xu, and J. Dong, "A review on supercooling of Phase Change Materials in thermal energy storage systems," *Renew. Sustain. Energy Rev.*, vol. 70, pp. 905-919, 2017.

[58] G. R. Solomon, S. Karthikeyan, and R. Velraj, "Sub cooling of PCM due to various effects during solidification in a vertical concentric tube thermal storage unit," *Appl. Therm. Eng.*, vol. 52, no. 2, pp. 505-511, 2013.

655

[59] Y. Korin, A. Fiterman, and G. Ziskind, "Modeling of Subcooling in Solidification Process."

[60] S. Höhlein, A. König-Haagen, and D. Brüggemann, "Thermophysical Characterization of MgCl2·6H2O, Xylitol and Erythritol as Phase Change Materials (PCM) for Latent Heat Thermal Energy Storage (LHTES)," *Materials*, vol. 10, no. 4, Apr. 2017.

660

[61] E. P. Del Barrio, R. Cadoret, J. Daranlot, and F. Achchaq, "New sugar alcohols mixtures for long-term thermal energy storage applications at temperatures between 70 C and 100 C," *Sol. Energy Mater. Sol. Cells*, vol. 155, pp. 454-468, 2016.

665
# Appendix

Detailed calculation of the matrices expression for the heat equation:

The equation (1) can be rewritten as a matrix system, considering the first order time derivation:

$$[G][T]^{k+1} + [Q] = \frac{1}{\Delta t} [diag(C_i)] ([T]^{k+1} - [T]^k) \quad (17)$$

It can be noted that in the boundary conditions terms here are not explicit. They can be expressed as van Neumann (as done here) or Dirichlet and are included in the source matrix [Q] and with the help of the conductance matrix via $G_{ii}$. The following equation is the development of the last expression.

$$([G] - \frac{1}{\Delta t} [diag(C_i)]) [T]^{k+1} + [Q] = - \frac{1}{\Delta t} [diag(C_i)] [T]^k \quad (18)$$

Then the isolation of the temperature vector at the next time step $k+1$ is done:

$$([G] - \frac{1}{\Delta t} [diag(C_i)]) [T]^{k+1} = - \left( \frac{1}{\Delta t} [diag(C_i)] [T]^k + [Q] \right) \quad (19)$$

This corresponds to the decomposed equations (2).

Detailed conditions for the supercooling factor:

Initialization

$$\begin{cases} \text{if } T(i, 0) < T_m + 2 \cdot \Delta T \quad \forall i & f_{\text{super}}(i, 0) = 0 \quad \forall i \\ \text{if } T(i, 0) \ge T_m + 2 \cdot \Delta T \quad \forall i & f_{\text{super}}(i, 0) = 1 \quad \forall i \end{cases} \quad \% \text{Initialization at } t = 0 \text{ for all PCM nodes } i \quad (20)$$

During the calculation loop

$$\begin{cases} \text{if } [T(i,t) \ge T_m + 2 \cdot \Delta T \quad \forall i] & f_{\text{super}}(i,t) = 1 \quad \forall i \\ \text{if } f_{\text{super}}(i,t) = 1 \quad \text{AND} \quad [T(i,t) \le T_c \text{ OR } f_{\text{super}}(\{i-1; i+1\}, t) = 0] & f_{\text{super}}(i,t) = 2 \\ \text{if } f_{\text{super}}(i,t) = 2 \quad \begin{cases} \text{if formulation } i = 1 & T(i,t) = T_{\text{eq}}^{(1)}(i) & f_{\text{super}}(i,t+\Delta t) = 0 \\ \text{if formulation } i = 2 & T(i,t) \ge T_{\text{eq}}^{(2)}(i), f_{\text{super}}(i,t+\Delta t) = 0 \end{cases} \end{cases} \quad (21)$$

Detailed calculation of the conservation of enthalpy

$$A_0 = A_1 + A_2 \Leftrightarrow A_0' = A_1'$$

$$C p^l (T_{eq}^{(1)} - T_c) = \int_{T_{eq}^{(1)}}^{+\infty} \frac{H}{\Delta T \cdot \sqrt{\pi}} \cdot \exp(-\frac{(T - T_m)^2}{\Delta T^2}) \cdot dT \quad (22)$$

$H$ is assumed to be constant so it can be removed from the integral expression.
$$Cp^l(T_{eq}^{(1)} - T_c) = \frac{H}{2} \int_{T_{eq}^{(1)}}^{+\infty} \frac{2}{\Delta T \cdot \sqrt{\pi}} \cdot \exp\left[-\left(\frac{T - T_m}{\Delta T}\right)^2\right] \cdot dT$$

$$Cp^l(T_{eq}^{(1)} - T_c) = \frac{H}{2} \cdot \frac{2}{\sqrt{\pi}} \int_{x_{eq}^{(1)}}^{+\infty} \frac{(T_{eq}^{(1)} - T_m)}{\Delta T} \exp(-x^2) \cdot dx$$

$$Cp^l(T_{eq}^{(1)} - T_c) = \frac{H}{2} \cdot \left( \frac{2}{\sqrt{\pi}} \int_{0}^{+\infty} \exp(-x^2) \cdot dx - \frac{2}{\sqrt{\pi}} \int_{0}^{x_{eq}^{(1)}} \exp(-x^2) \cdot dx \right)$$

$$Cp^l(T_{eq}^{(1)} - T_c) = \frac{H}{2} \cdot \left( 1 - \frac{2}{\sqrt{\pi}} \int_{0}^{x_{eq}^{(1)}} \exp(-x^2) \cdot dx \right)$$

$$Cp^l(T_m - T_c) + Cp^l(T_{eq}^{(1)} - T_m) = \frac{H}{2} \cdot \left( 1 - \operatorname{erf}\left(\frac{(T_{eq}^{(1)} - T_m)}{\Delta T}\right) \right)$$

685 Assuming $\Delta T$ is small enough so $Cp^l(T_{eq} - T_f)$ is negligible.

$$Cp^l(T_m - T_c) = \frac{H}{2} \cdot \left( 1 - \operatorname{erf}\left(\frac{(T_{eq}^{(1)} - T_m)}{\Delta T}\right) \right)$$
$$\operatorname{erf}\left(\frac{(T_{eq}^{(1)} - T_m)}{\Delta T}\right) = 1 - \frac{2 \cdot Cp^l(T_m - T_c)}{H} \qquad (23)$$

$$T_{eq}^{(1)} = T_m + \Delta T \cdot \operatorname{erf}^{-1}\left(1 - \frac{2 \cdot Cp^l(T_m - T_c)}{H}\right)$$

<table><caption>Crystallization rate conditions for the different cases in Figure 12</caption><thead><tr><th></th><th>a<sub>0</sub> (-)</th><th>a<sub>1</sub> (K<sup>-1</sup>)</th><th>a<sub>2</sub> (K<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Configuration #1</td><td>0.12</td><td>0</td><td>0</td></tr><tr><td>Configuration #2</td><td>0.012</td><td>0</td><td>0</td></tr><tr><td>Configuration #3</td><td>0.007</td><td>0.003</td><td>0</td></tr><tr><td>Configuration #4</td><td>-0.055</td><td>0.003</td><td>0</td></tr></tbody></table>

<caption>Table 3: crystallization rate polynomial coefficients for the test of different characteristic laws</caption>