One-Dimensional Systemic Modeling of Thermal Sensors Based on Miniature Bead-Type Thermistors

Accurate measurements of thermal properties is a major concern, for both scientists and the industry. The complexity and diversity of current and future demands (biomedical applications, HVAC, smart buildings, climate change adapted cities, etc.) require making the thermal characterization methods used in laboratory more accessible and portable, by miniaturizing, automating, and connecting them. Designing new materials with innovative thermal properties or studying the thermal properties of biological tissues often require the use of miniaturized and non-invasive sensors, capable of accurately measuring the thermal properties of small quantities of materials. In this context, miniature electro-thermal resistive sensors are particularly well suited, in both material science and biomedical instrumentation, both in vitro and in vivo. This paper presents a one-dimensional (1D) electro-thermal systemic modeling of miniature thermistor bead-type sensors. A Godunov-SPICE discretization scheme is introduced, which allows for very efficient modeling of the entire system (control and signal processing circuits, sensors, and materials to be characterized) in a single workspace. The present modeling is applied to the thermal characterization of different biocompatible liquids (glycerol, water, and glycerol–water mixtures) using a miniature bead-type thermistor. The numerical results are in very good agreement with the experimental ones, demonstrating the relevance of the present modeling. A new quasi-absolute thermal characterization method is then reported and discussed. The multi-physics modeling described in this paper could in the future greatly contribute to the development of new portable instrumental approaches.


Introduction
The miniaturization of electro-thermal sensors consists of designing and manufacturing ever smaller measuring devices, in order to characterize very limited quantities of material, under a wide variety of operating conditions. This miniaturization has opened up new perspectives in many scientific and technical fields, particularly in medicine, where it has enabled in vivo thermal characterization of biological tissues, that is of crucial importance, for example, in the focal treatment of tumors by magnetic hyperthermia [1].
If some attempts of thermal systems miniaturization at the microscopic scale have been successfully achieved, particularly with the development of labs on a chip, they mainly concern the detection of micro objects [2], thermal imaging, or heat transfer enhancement in micro-fluidics [3]. The precise measurement of the thermal properties of materials is still mainly carried out today at centimeter or millimeter scales [4,5]. Therefore, in this paper, we refer to miniature sensors as sensing elements whose main dimensions are in the order of the millimeter to a few hundred micrometers.
It is clear that the miniaturization of electro-thermal sensors brings undeniable advantages, such as: a lower thermal inertia; a limited invasiveness of the system to be characterized, which is, therefore, less disturbed by the sensor; and the possibility of using small quantities of material, that is highly desirable in nanotechnologies, for example, where the products to be characterized can be very expensive. However, the miniaturization of electro-thermal sensors can also have drawbacks, the most delicate being certainly the problem of the increasing influence of boundary-effects, mainly due to electrical connections. These connections represent unavoidable heat sinks for the sensor, systematically imposing a precise and delicate calibration before any use [6,7]. It is mainly this aspect that makes it very difficult to design new measuring systems using miniature electro-thermal sensors. Under these conditions, it is very desirable to be able to develop systemic modelings, allowing for describing, as accurately as possible, the real operation of thermal sensors in their measurement environment. Thus, a one-dimensional (1D) systemic modeling using a Godunov-SPICE discretization is presented in this paper, in the case of miniature bead -type thermistor sensors with negative temperature coefficient (NTC). This approach allows for a very efficient model of the sensor, its control, and signal processing circuits, as well as the thermal interactions with the material to be characterized, and this within a single workspace (systemic modeling).
There are a large number of thermal characterization methods (depending on the nature of the materials to be characterized), which can be classified mainly into two categories [8][9][10], without being exhaustive: • the transient methods, among which we can mention: the pulse-decay method [11]; the temperature step methods [6,[12][13][14]; the transient hot wire (THW) and the transient hot strip (THS) methods; flash methods; transient plane source methods (TPS); 3ω harmonic methods [15,16]; • the steady-state methods, such as the guarded hot plate and the heated thermocouple methods.
These techniques can be implemented using a wide variety of different thermal sensors (temperature or heat flow), such as: thermocouples, negative temperature coefficient (NTC) thermistors, resistance temperature detectors (RTDs), and radiative sensors.
In the present paper, we will focus on the transient mode use of miniature resistive electro-thermal sensors, for thermal characterization of solid or liquid materials, in the laboratory or in everyday use conditions. These sensors can be classified into two main categories: the hot metal wires (HMW) or hot metal films (HMF) made of pure metal resistors (RTD), such as platinum or nickel; and the negative temperature coefficient (NTC) thermistors, generally composed of semiconducting metal oxides. All these sensors share the same operating principle: the time evolution of their electrical resistance R(t), due to a controlled self-heating, depends on the physical properties (thermal conductivity, density, specific heat, fluid flow) of the material to be characterized, in which the sensor is immersed.
Among these resistive electro-thermal sensors, bead-type thermistors have the advantages of spherical or quasi-spherical symmetry, being non-traversing and small in size, unlike wire sensors (which are traversing or large in size) and hot film sensors (which are usually non-symmetrical composite structures [17]). In contrast, RTDs have the great advantage of obeying a law of variation of the resistance as a function of temperature which is linear with a very good approximation, while NTCs obey a non-linear law of variation. The linearity property of RTDs greatly simplifies the development of signal processing circuits and allows a wide variety of applications, ranging from thermal characterization of materials (THW method, 3ω method) to flow measurement (hot wire or hot film anemometry).
In contrast, the non-linearity of NTCs has certainly been a major obstacle to the development of various thermal characterization methods using them. If NTCs are indeed used for transient thermal characterization of pulse-decay or temperature step types, however, we do not find in the literature any use in harmonic method of 3ω type for example. Thus, main objective of this paper is to propose an accurate 1D systemic multi-physics modeling of bead-type thermistors, which are non-linear but very easy to implement sensors. This systemic modeling should, on the one hand, help in the development of new thermal characterization methods using miniature bead-type NTCs; on the other hand, it should help to better understand the operating mode of existing methods.
The implementation modes of miniature resistive sensors used in material characterization can be classified into two main categories: • the source/sensor(s) mode [18][19][20], which implements separate heat source and sensors. This operating mode requires an accurate and reproducible positioning of the sensor(s) relative to the heat source; and • the self-heating mode [10], in which the sensitive element serves as both a heat source and a temperature sensor. The basic principle of this mode of operation is that the less heat conducting the surrounding materials to be characterized, the greater the self-heating of the sensor element. Thus, the average temperature T c (t) of the sensor and the electrical power P e (t) it dissipates are, in most situations, the basic informative signals considered in electro-thermal methods used for thermal characterization.
While the 1D systemic modeling proposed in this work allows the simulation of both characterization modes, it is the self-heating mode that will be mainly exposed here.
Simulation has become an essential tool in the development of electronic circuits. In electronics, the SPICE (Simulation Program with Integrated Circuit Emphasis) code, developed at Berkeley in the early 1970s [21], has been widely adopted and has indeed become a true industrial standard. While SPICE excels at solving the ordinary differential equations (ODEs) that describe the temporal operation of electrical circuits, it lacks the intrinsic ability to solve the partial differential equations (PDEs) that describe heat transfers and spatiotemporal evolutions of the temperature within the sensor and the material being characterized. To correct this deficiency, it is necessary to develop accurate models of electro-thermal sensors, which integrate interactions between thermal quantities and electrical signals. In reference [22], Heyd et al. used the capability of SPICE for fast and accurate Laplace transform inversion to solve the unsteady heat equation in a hot wire heated by Joule effect and submitted to a fluid flow (hot wire anemometer). In this context, time and position were treated as continuous variables, allowing almost unlimited accuracy in the simulation of the different modes of operation of hot-wire anemometers. Other approaches are possible [23][24][25][26][27][28][29], which generally implement, on the one hand, a spatial discretization of the medium to be characterized (and possibly also of the sensor) and uses, on the other hand, an analogy to transform this discretization into an analogous electrical circuit. This analogous circuit is then coupled to the electrical control signals.
A hybrid approach is used in the present work. The material to be characterized and the composite bead-type NTC (core and sheath) are all discretized in spherical symmetry (1D equivalent model) using a Godunov-type numerical scheme [30], while the electrical connection wires of the thermistor are modeled using thermal resistors. The results obtained using this numerical model in SPICE are then compared with the experimental results obtained with pure glycerol at rest. Then, different effective parameters of the 1D equivalent model are adjusted to allow the best possible reproduction of the experimental measurements.

Glycerol
In order to avoid natural convection as much as possible during the development of the 1D systemic NTC model, glycerol, a non-ionic liquid with high viscosity, was chosen as test fluid.
Because of the high hygroscopicity of glycerol, only fresh liquid (Sigma-Aldrich, St, Saint-Quentin-Fallavier, France) was used, and sealing precautions were taken to avoid dissolution of atmospheric water vapor into the liquid during measurement and storage.

NTC Bead-Type Thermistor
Accurate and reproducible determination of the thermal properties of materials using a miniature NTC bead-type thermistor requires a good understanding of the different heat exchanges that take place within the system. It is difficult, if not impossible, to consider absolute measurements of thermal conductivity and diffusivity using NTC thermistors, without a precise knowledge of their internal structure and of the influence of electrical contacts.

Composition-Electrical Properties
The sensor used in this work is a miniature high precision epoxy encapsulated beadtype thermistor of type 44004 from Omega, whose resistance at a temperature of 25 • C is equal to 2252 Ohms (2252@25; see Figure 1a). This type of high precision thermistor was chosen both for its quasi-spherical shape and small size and, on the other hand, for its interchangeability. To avoid electrical conduction problems in ionic solvents, such as water, a thin layer of insulating varnish has been applied to the NTCs copper connecting wires. Miniature NTC thermistors of bead-type generally consist of an active core (see Figure 1), composed of a semiconducting metal oxide with negative temperature coefficient, protected by a glass, nylon or epoxy sheath. While the active core can sometimes have a spherical shape, it is, rather, usually in the shape of a slab (see Figure 1b). In order to allow the circulation of an i(t) electric excitation current through the core, two small metal electrodes are deposited on it, also allowing for measurement of its electric resistance R(T), which depends on the core average temperature T.
The external shape of the NTC is usually a prolate spheroid, whose characteristic dimension d e is of the order of the millimeter or less (2.4 mm in the case of 44004 model). Two small diameter copper wires (in the order of a few hundred micrometers: 180 µm for 44004) ensure the electrical connection of the active core with the control and processing circuit of the electro-thermal signal.
To measure the temperature T of the NTC active core, it is necessary to be able to relate T to the value of the NTC electrical resistance R(T). There are several common laws used to describe the variations in the electrical resistance R (in Ω) of the active core as a function of its equilibrium temperature T (in K). One of the most used is certainly the law of Steinhart and Hart [31]: where the constants c 0 , c 1 , and c 3 are obtained by a precise calibration (at least three points) of the thermistor considered, in the temperature range visited experimentally. In the case of thermistor 2252@25, for example, the values of R(T) (taken from the "Thermistor Resistance vs. Temperature" table provided by Omega) in the range [−80, 120] • C allow the following coefficients to be found (in K −1 ): c 0 = 1.47 × 10 −3 , c 1 = 2.38 × 10 −4 , and Simpler models can also be used [32], with only two coefficients. For example: where T ref is typically chosen equal to 298.15 K. In the case of the 2252@25 thermistor, we have obtained R ref = R(T ref ) = 2252 Ω and β = 3864.5 K. This kind of two-factor relationship being less precise than (1); therefore, it is recommended to experimentally determine the coefficients (b 0 , b 1 ) or (R ref , β) in the range of temperatures investigated during the thermal characterization. We now review the main energy transfers that take place within the thermistor and its surrounding environment. During a typical liquid thermal characterization experiment, the average electrical power P e dissipated by the thermistor is in the range of 1 to 10 mW. The corresponding temperature variations ∆T involved are then in the order of a few Kelvins, thus generally being much lower than the working temperatures.
Considering the 2252@25 Omega thermistor used in this work, having an average external radius a e ≈ 1.2 mm, typically traversed during ∆t = 30 s by an electric current I e = 2 mA, this corresponds to an average dissipated electric power P e ≈ 10 mW. When this thermistor is immersed in pure glycerol (at rest), at the working temperature T 0 = 294 K, the maximum temperature variations in the active core are about ∆T c = T c − T 0 ≈ 1.3 K and ∆T e = T e − T 0 ≈ 0.8 K for the outer surface of the protective sheath. Considering these orders of magnitude typical of a thermal characterization experiment, we can now analyze the main sources of thermal losses and poor thermal contacts: • Radiative heat losses P r at the external surface of the NTC: These losses can be estimated using Stefan-Boltzmann's law: . Under the working conditions presented above, and considering a maximal emissivity ε = 1, we obtain P r ≈ 0.1 mW. Therefore, it is justified to neglect P r in front of P e . Under the typical working conditions described in this work, it can be assumed that radiation losses have a negligible influence on the thermal characterization of materials using a miniature bead-type thermistor. • Thermal losses P L at the electrical contacts between the thermistor and the control circuit: It is difficult to precisely estimate these heat losses because they are closely linked to the nature of the experimental system implemented: electrical insulation of the connecting wires (using a varnish or silicone) and length L of the connecting wires; use of metal or plastic mounting probe. Under the operating conditions described above, and for varnish insulated copper wires, with a diameter d w = 0.18 mm and a length of L = 1 cm, immersed in pure glycerol at rest, a finite element modeling has led to heat losses typically in the order of 3 mW. Therefore, heat losses via the electrical connections are considerable and can in no way be neglected in the modeling of the thermistor. • Convective losses P c when the material to characterize is a fluid: Convective losses always lead to an additional heat extraction from the sensor and are at the origin of an overestimation in the measurement of the conductivity of the fluid to be characterized. Therefore, it is essential to limit these convective losses if we want to correctly estimate the thermal conductivity of the fluid to be characterized. This requirement dictated the choice of glycerol as a test liquid for 1D systemic modeling of NTCs. The thermal power evacuated from the thermistor to the fluid can be evaluated in the presence of natural convection, from the Churchill correlation [33], which gives the following Nusselt number expression, valid for natural convection around a sphere: where Pr = ν/α > 0.5 is the Prandtl number of the fluid, and Ra d e = gβd 3 e ∆T e /να < 10 11 is the Rayleigh number for the bead-type thermistor, with: g = 9.81 m · s −2 the acceleration of gravity; ν = η/ρ the kinematic viscosity of the fluid (m 2 /s); α = k/ρc its thermal diffusivity (m 2 /s); and β its volumetric expansion coefficient (K −1 ). Note that, in the absence of natural convection (Ra d e = 0), the Nusselt number takes the limiting value Nu = 2, which corresponds to pure thermal conduction through the liquid at rest. Typical values (at 20 • C) of water physical properties are collected in Table 1. Typical values of glycerol properties at 20 • C can be found for example in [34]. Using the values of Table 1 for ∆T e = 0.8 K, we find in the case of glycerol that Nu ≈ 2.5, which remains close to the value obtained in the absence of natural convection, and P c = 0.8 mW, which is negligible compared to P e . It can be concluded that it is acceptable, for temperature variations ∆T e 1.0 K near room temperature, to neglect the contribution of natural convection around the NTC in the case of immersion in glycerol. In contrast, in the case of NTC immersion in water, considering a temperature difference ∆T e = 0.5 K (thermistor heating is less here than in the case of glycerol) near room temperature, we find that Nu ≈ 3.7, which is quite different from the values obtained in the absence of natural convection, and P c = 3.8 mW, which can no longer be ignored. • Poor thermal contacts: There are two main sources of poor thermal contact here: the thermal contacts between the active core of the NTC and its protective sheath and between the sheath and the surrounding medium to be characterized. If the latter is a fluid (which is the case in the present study), we can assume that the corresponding sheath/fluid thermal contact resistance is negligible. In contrast, the contact resistance R c between the sheath and the core is potentially significant, unknown and a priori different from one sensor to another and will be taken into account in the 1D systemic modeling proposed in this study. Despite its small size, a miniature bead-type NTC is a complicated object from a heat transfer point of view. Indeed, the presence of both metallic electrodes and connecting wires, but also of an active core that can be slab-shaped and surrounded by a prolate spheroid protective shell, does not strictly allow for attribution of particular simplifying symmetries to the T(M, t) system temperature field. As a result, an exact description of the thermal behavior of the "Fluid+ NTC + control circuits" multi-physics system would require a three-dimensional (3D) modeling of heat transfers within the system. Unfortunately, the complexity of such a 3D finite element modeling is a significant obstacle to the development of an efficient and fast systemic modeling of NTCs used in thermal characterization.
However, and given the small size of the miniature bead-type thermistor relative to the volume of fluid to be characterized, the heat transfers outside the active core depend very little on its shape (sphere or slab), as illustrated by the finite element modeling results presented in Figure 2. These results were obtained by using the same constant and uniform heating power densityq e in both cases. The effective radius a c of the equivalent spherical active core (see Figure 2b) was adjusted here to obtain the same core temperature as in the slab-shaped core case (see Figure 2a). Note that, in this case, the volume V c = 4 3 πa 3 c of the equivalent sphere is equal to the volume of the slab.
(b) Spherical active core. Finite elements modeling of heat transfers through the multi-physics system: NTC + Fluid, with: (a) a slab-shaped active-core; (b) a spherical active core. In both cases, the same materials were considered for the sheath (epoxy, external radius a s = 1 mm), the active core (semiconductor material), and for the surrounding fluid (radius a f = 5 mm) to be characterized.
Based on the results presented in Figure 2, it can be concluded that, regardless of the shape (sphere or slab) of the NTC active core, it is possible to propose an approximate 1D modeling of a miniature bead-type NTC, which uses an equivalent spherical active core of effective radius a c . This simplification allows us to propose a quasi-1D modeling of the NTC, with spherical symmetry. Note that, in the case of spherical core NTCs, the effective radius a c can be identified with the exact radius of the active core.

Mathematical Modeling
An elementary mathematical model, that is frequently used to describe the thermal behavior of a miniature NTC of bead-type, consists of assimilating the whole composite sensor to a single sphere of temperature T c (r, t), with an effective radius a, and constant effective thermal properties: k (thermal conductivity) and α (thermal diffusivity) [12]. The self-heating of the sensor by the Joule effect is usually modeled by a uniform power densitẏ q e = P e (t)/ 4 3 πa 3 c , where P e (t) is the electrical power supplied by the driving circuit to the active core of the NTC. The heat transfers through the NTC thermistor, being solely of a diffusive nature, and its temperature T c (r, t), thus, obeys the following heat equation in spherical symmetry [36], with r being the distance to the center of the sphere: with T c (0, t) finite and T c (r, 0) = T 0 . According to the results and discussions of Section 2.2.2, it is also considered in the present work that the thermistor is spherical in shape, of radius a s , but it is decomposed into two concentric parts (see Figure 3): the active semiconducting core (Figure 3c) of effective radius a c < a s , thermal conductivity k c , and temperature T c (r, t); the insulating protective sheath (s) of effective thickness e s = a s − a c , thermal conductivity k s k c , density ρ s , specific heat c s , and temperature T s (r, t). As in the case of the model proposed by Balasubramaniam et al. [12], the present ideal model does not take into account at this stage the losses due to the connecting wires. These defects will be modeled later using appropriate thermal resistances.
r active core Figure 3. Ideal 1D effective model of a miniature bead-type NTC, immersed in a fluid at rest (with thermal conductivity k f , density ρ f , and specific heat c f ). The red curve represents effective variations of the δT = T(r, t) − T 0 temperature through the system, as a function of r, at a given time t.
Due to the usual high value of the thermal contrast (or thermal ratio k c /k s ) between the core and the sheath of the NTC, the temperature T c (r, t) = T 0 + δT c (r, t) of the high thermal conductivity core can be considered as uniform to a very good approximation (see Figure 2). Numerical simulations made with the multi-physics finite element solver Flex-PDE have shown that, in all the cases under study here, the spatial variations in the temperature δT c (r, t) of the core were absolutely negligible: Thus, it is possible to consider here that the temperature δT c only depends on the time t. Therefore, the energy balance of the highly conductive core can be written, using the spherical symmetry of the equivalent 1D model, as: where ρ c and c c are the density and the specific heat of the active core, respectively, and T s is the temperature of the sheath, at the interface with the active core at r = a c . On the other hand, the temperature T s (r, t) of the sheath obeys a diffusive heat equation without source term, which is written, still using the spherical symmetry of the equivalent 1D model, as: where the thermal diffusivity α s of the sheath is supposed to be constant. The thermal contact between the solid sheath (s) and the solid core (c) being a priori not perfect, we can only assume here the continuity of the heat flow at the interface r = a c between (c) and (s): The surrounding material (f) to be characterized is assumed to be homogeneous, at rest, of semi-infinite extent and constant thermal conductivity k f and diffusivity α f . Two situations are often considered, depending on the physical nature of the surrounding material, delimited by the radius r = a f : is an inert (non-biological) material at rest, with no heat source term. In this case, its temperature T f (r, t) also obeys a diffuse heat transfer equation, with spherical symmetry: • (f) is a biological material, which obeys a differential equation of the Penne type (bioheat transfer equation), for example [37,38]: whereq f and ωρ b c b (T 0 − T f ) are the metabolic and the perfusion heat source terms, respectively. The perfusion term describes the heat transfer between the medium (f) (biological tissues) and blood (b) flowing through the veins and arteries, with a volumetric mass flow ρ b ω (kg·m −3 ·s −1 ) and blood specific heat c b . The initial condition to be considered for the partial differential Equation (8) or (9) is The boundary condition to consider at r = a f depends on the thermal characterization apparatus that is used. We suppose here that the medium to characterize is in perfect thermal contact at r = a f with a thermostat at the temperature T 0 : While both situations (8) and (9) can be studied in the same way with the systemic modeling presented in this work, it is the (8) case that will be presented here in detail, both from a numerical and an experimental point of view. Finally, a perfect thermal contact is assumed between the thermistor sheath and the medium to be characterized (fluid at r = a s ), the following continuity relationships must then be satisfied at all times: The simultaneous solving of the (PDEs) partial differential Equations (5), (6) and (8) (or (9)), together with the boundary conditions (7), (10), (11), and the initial condition T c (0) = T s (0, t) = T f (0, t) = T 0 , allows for determination of the temporal evolution of the core temperature T c (t), which depends both on the medium to be characterized and the electrical power P e (t) used to excite the thermistor.
Thus, the 1D systemic approach presented in this work consists of solving the PDEs of the equivalent system, while, at the same time, taking into account the operation of the electronic circuits used to provide the power P e (t) to the thermistor.
The electronic circuit that was used for both the experimental measurements and the systemic modeling is described in detail in the following paragraph.

Electronic Circuit
The electronic circuit shown in Figure 4 was used for both experimental measurements and systemic modeling. This is a typical example of an amplified voltage divider bridge circuit, quite versatile, which can be adapted according to the implemented resistive sensor (NTC or RTD) and can be used both in transient and frequency modes.
A Measurement Computing (MC) USB-2537 multi-functions data acquisition board has been used for measurements, providing one 16 bits DAC output and two 16 bits ADC inputs (ADC0 and ADC1). The following precautions must generally be observed when implementing transient thermal characterization methods with the circuit shown in Figure 4: • Since the temperature variations due to self-heating are generally of small amplitude (a few Kelvins at most), it is necessary to use an analog-to-digital converter (ADC) with a resolution of at least 14 bits. Moreover, as thermal phenomena are generally relatively slow (with time constants typically in the order of a few tenths of a second or more), a sampling period of around 10 ms or more is, therefore, sufficient in most cases. A digital-to-analog converter (DAC) with a resolution of 12 bits or better can be used to provide the excitation signal v e (t) for transient methods. If the sensing element is to be excited in current rather than voltage, a voltage-to-current converter (of the Howland source type) can be used instead of the non-inverting amplifier shown in Figure 4. • The inputs ADC0 and ADC1 must have a sufficiently large input impedance Z e (typically, Z e > 10 5 Ω. In the case of USB-2537: Z e = 10 7 Ω) to not load the circuit. If necessary, a high impedance buffer (voltage follower) can be used to isolate the control circuit from the influence of the ADC. • The operational amplifiers (OA) used in circuit Figure 4 must operate in the linear regime. If the self-heating of the sensor requires an electric current with an intensity i(t) > 10 mA, then, the use of an operational amplifier capable of delivering high currents should be considered. This could be the case, for example, with low-resistance platinum wires, whose value is close to 1 Ω. In this case, a typical average electrical power P e = 10 mW requires about 100 mA electrical excitation current.
To make the set-up as versatile as possible, a power OA of type L272 (delivering currents up to 1 A without significant harmonic distortion) was systematically used. • The working (or baseline) temperature T 0 must be precisely regulated, usually by means of a temperature controlled bath. The variations ∆T 0 of the working temperature must be negligible in front of the maximal temporal variations of the sensor core temperature δT c (t), due to self-heating. An accuracy of ∆T 0 ≈ 0.01 K is sufficient in most cases. Figure 4. Amplified divider bridge circuit used both to control the self-heating of the NTC bead-type thermistor and measure the time variations of its electrical resistance R(T c ). This circuit can be used in both transient and frequency modes. ADC0 and ADC1 are two inputs of an Analog to Digital Converter (ADC). The circuit can be excited either by using a Digital To Analog Converter (DAC) or by using an Arbitrary Waveforms Generator (AWG).

1D Systemic Modeling
As mentioned before, an approximate 1D systemic modeling of a bead-type NTC is considered in this work, which consists essentially of assuming that heat transfers through the NTC present the spherical symmetry. This 1D systemic modeling of a bead-type NTC relies first on a suitable spatial discretization of the heat Equation (6) in the sheath and (8) or (9) in the medium to be characterized.
The discretization process of the present thermal problem, as well as different elementary applications of this discretization, is introduced in the present paragraph.

General Approach
In all cases considered in the present study, we are dealing with (equivalent) domains having the shapes of spherical shells, for which the discretization schemes are, therefore, perfectly similar. Therefore, we choose to discretize the partial differential Equation (9), which has the most general expression here. This PDE is then rewritten in the following general self-explanatory form: where α = k/ρc is the thermal diffusivity of the physical domain (Ω) = (a 1 , a 2 ). The PDE (12) can be rewritten in the following form: where Figure 5).
Using explicit discretization schemes centered on the node r m = a 1 + mδ r , the discrete version of (13) can be written as follows: where T n m is the absolute temperature of the node (r m ) at time t n = t n−1 + ∆t n with n 1; m ∈ [1, M − 1]. Equation (14) can now be rewritten as a thermal energy conservation relation (Godunov conservative form), where the electro-thermal analogies gathered in Table 2 were used: i th

Thermal Electro-Thermal
Thus, it is possible to give an electro-thermal representation of the slice (m), which is composed of five ideal controlled current sources and one capacitor (see Figure 6). Note that the additional ideal voltage source guarantees an initial temperature of node (r m ) equal to T 0 . Moreover, the capacitance C th of the analog electro-thermal capacitor is given by C th = ρcδ r .
h ω m slice (m) Figure 6. Electro-thermal modeling of a slice (m) of a material obeying the partial differential Equation (12).
To test the relevance of the Godunov electro-thermal analogy proposed in this work, we first solved Equation (12) in the case whereq = 0 and h ω = 0, considering the following initial and boundary conditions: δT(r, 0) = T(r, 0) − T 0 = 0, δT 2 = T(a 2 , t) − T 0 = 0, and δT 1 = T(a 1 , t) − T 0 = 10 K (Dirichlet conditions). The considered system is a spherical shell of inner radius a 1 = 10 −3 m and outer radius a 2 = 2 × 10 −3 m, with the following physical properties: k = 1.0 W · m −1 · K −1 ; ρ = 1.9 × 10 3 kg · m −3 and c = 1.5 × 10 3 J · kg −1 · K −1 . The spherical shell is discretized into M = 12 slices of equal thickness (see Figure 7a). Figure 7b shows the results obtained with the present systemic modeling (symbols) when t → ∞ (steady state), for which case there is an analytical solution, given by: It can be seen from Figure 7b that the numerical results are in very good agreement with the analytical results (solid line). Figure 7c shows the variations of temperature δT(r m , t) as a function of time at different points r m throughout the shell. The results obtained with the present modeling match remarkably well with those provided by a numerical resolution using the finite elements method (FlexPDE 6.0).
The case of a constant heat volumetric sourceq = cste was also computed, with the other properties being identical to those used previously. There is still an analytical steady-state solution in this case, and, considering the same Dirichlet boundary conditions as used previously, one gets: Figure 8a (steady state profile) and Figure 8b (temporal variations at different points) show the results obtained with the present approach forq = 10 8 W · m −3 and M = 12. Figure 8a compares the results obtained with the present electro-thermal analogy to the analytical solution (17) for the case of steady state profile, while Figure 8b compares the results of the electro-thermal analogy to those obtained by a finite elements solver. In both cases, it can be seen that the correspondence between the two approaches is excellent. Using current sources at the boundaries rather than voltage sources (see Figure 9a), it is also possible to simulate Neumann-type boundary conditions: Φ(a i ) = (∂ r T) r=a i . In the case of a Dirichlet boundary condition at r = a 1 and a Neumann boundary condition at r = a 2 , there is an analytical expression of the steady state solution, which is given here by: The numerical resolution being more demanding with this kind of boundary condition, it is generally preferable to use a more precise discretization in this case. Figure 9b (steady state profile) and Figure 9c (temporal variations at r = a 2 ) compare the results obtained forq = Φ(a 2 ) = 0 using the present approach with M = 32, to those obtained by a finite elements solver, with the other properties being identical to those used previously. In both cases, it can be seen that the correspondence between the two approaches is still very good.  Finally, solving the bioheat Equation (13) is also considered, with T 0 = 0 K and using Dirichlet boundary conditions: T(a 2 , t) = T 2 = 0 K and T(a 1 , t) = T 1 = 1 K. Usual values of the physical parameters were chosen from Bergman et al. [38] for human tissues and blood:q = 700 W · m −3 , h ω = 1800 W · m −3 · K −1 , ρ = 1000 kg/m 3 , c = 3600 J · kg −1 · K −1 , and k = 0.5 W · K −1 · m −1 . Note that an analytical expression also exists here for the steady-state temperature, which is given by: where β = √ h ω /k. Figure 10a compares the results obtained by the present electro-thermal analogy with M = 28, to the analytical solution (19) for the case of steady-state profile. The maximum relative error between numerical and analytical solutions is in this case less than 6 × 10 −4 . On the other hand, Figure 10b compares the variations of T(a 1 + δ r , t) as a function of time, given by the electro-thermal analogy, to the variations obtained by a finite elements solver. In both cases, it can be seen that the correspondence between the two approaches is still excellent.

NTC Bead-Type 1D Systemic Modeling
The previous examples have shown that the approach developed in this work can handle a wide variety of spherically symmetric heat transfer problems, close to the problem of thermal characterization of materials using a bead-type NTC. Thus, we are now ready to propose a 1D systemic modeling of a bead-type NTC, immersed in a fluid at rest.

Ideal NTC Bead-Type 1D Systemic Modeling
We first consider the ideal case for which the thermal contact between the active core and the sheath is perfect, and the thermal losses along the connecting copper wires are not taken into account. In order to be consistent with the electro-thermal current densities gathered in Table 2, the energy balance (5) of the active core is divided by the effective surface area of the equivalent spherical core S c = 4πa 2 c : In the case of an ideal thermal contact between the active core and the sheath, the second term of the right-hand side of (20) is already taken into account in the discretization of the sheath (see Figure 6, node m = 1). In this case, T c (t n ) ≡ T n 0 and the electro-thermal circuit shown in Figure 11 can be used to describe the temporal evolution of the ideal active core uniform temperature T c (t). wires sheath ideal active core model 1D equivalent core model Figure 11. One-dimensional electro-thermal systemic modeling of an ideal NTC active core with electrical resistance R(T c ) given by Equation (1) or (2). Recall that the effective radius a c is chosen so that the volume V c of the equivalent spherical core is equal to that of the slab.
We have found that simulating thermal characterization of materials by self-heating mode, using the ideal NTC core equivalent model shown in Figure 11 led to a systematic and significant overestimation of the thermistor active core temperature rise. This overestimation is due to the fact that the ideal model does not take into account thermal losses along the copper wires which connect the active core at temperature T c to the measurement and control circuit at working temperature T 0 < T c . In practice, these thermal losses limit considerably the temperature rise of the thermistor core, so they must be taken into account if an acceptable systemic modeling of the NTC is to be obtained.

Realistic NTC Bead-Type 1D Systemic Modeling
The ideal model shown in Figure 11 has been corrected by introducing an additional thermal path going from the thermistor core to the measurement circuit. This thermal path has been modeled here by the mean of a thermal resistance R L , whose value has been adjusted using experimental results obtained from the thermal characterization of pure liquid glycerol at rest (see Section 4.2). The NTC active core and its insulating shell being both solid bodies, the thermal contact between these two components cannot be strictly ideal. So, a thermal contact resistance R c between the core and the sheath has, therefore, also been introduced into the model. The value of R c , which was found to be much less critical than the value of R L , was also deduced from the thermal characterization of pure glycerol.
Finally, the complete 1D systemic model proposed in this work is represented in Figure 12. It takes into account the following physical properties and phenomena: an uniform temperature of the active core, heat transfers between the core and the sheath and between the sheath and the fluid to be characterized, the spatial discretization of the insulating sheath and of the fluid (in spherical symmetry), thermal paths due to the connecting wires, the contact thermal resistance between the sheath and the active core, and interactions between the NTC and electrical control circuits (through the electrical current i(t) and the electrical resistance R(T c )).
In order to adequately model the experimental conditions, Dirichlet boundary condition (BC) T f = T 0 was imposed on the fluid at r = a f . The sheath and the fluid elements of Figure 12 were both discretized using the electro-thermal model of Figure 6 and the electro-thermal current densities of Table 2. Note that it would also be possible to propose a more elementary 1D model than the one developed here. In this more basic model, the active core would be assumed to be punctual, and the heat transfer through the sheath would be modeled using a simple thermal resistance R s . The other elements of this point model would be similar to those used in the present 1D model.

Results
In this last part, experimental results obtained with different 44004 type NTCs are compared to the numerical results provided by the 1D systemic modeling developed in the previous paragraphs.
Glycerol has been chosen here as a reference liquid in order to determine the values of the different effective parameters needed by the 1D NTC spherical equivalent model.

Presentation
As mentioned in Section 2.1, high dynamic viscosity (η f = 1.48 Pa · s at 20 • C [34]) and thermal conductivity close to that of most biological fluids (k f = 0.285 W · m −1 · K −1 at 20 • C [34,39]) made liquid glycerol a good candidate to set the values of the physical parameters entering the 1D systemic modeling of a self-heated NTC used in thermal characterization mode. We recall that a thin layer of insulating varnish has been applied to the NTCs copper connecting wires.
Since NTC manufacturers give very little information about the physical properties of the materials used to produce their thermistors, typical starting values were assigned to a number of unknown parameters used in the model. These values were further refined by minimizing the following quantity, over a large number of tests (N test 20): where: T com c,n and T mea c,n are the computed and measured active core temperatures, respectively, at time t n = t n−1 + δt n ; N f is the number of experimental samples recorded during one test.
The model parameters values leading to the minimization of δT rms c for the entire set of tests performed are gathered in Table 3. Note that: • the model does not require the knowledge of the active core thermal conductivity k c because the core temperature T c is assumed to be uniform (k c k s ); • the values of k s and ρ s c s shown in Table 3 are compatible with usual epoxy values; and • the value of the shell external radius a s was determined using a caliper. The NTC being a prolate spheroid, an average value was considered.

Glycerol Thermal Characterization
A number of tests have been carried out with glycerol, which have led to the values shown in Table 3. The tests presented in this paper concern the immersion length of the NTC copper connecting wires.

Immersion Length Tests
In order to test the relevance of the loss resistance R L , which is a key parameter of the present systemic model, we have studied the temporal variations of the active core temperature T c (t) considering four different immersion lengths L i of the same thermistor in the same liquid (glycerol) at rest: L i = 0.0, 0.1, 5.0 and 10.0 mm (see Figure 13). The immersion of the NTC in the liquid was carried out using the mobile arm of a digital Vernier caliper providing a precision to 0.01 mm. A miniature digital microscope (with magnification ×10 to ×300) was used to accurately determine the position of the top of the NTC relative to the free surface of the liquid being characterized. The control circuit is that of Figure 4 Three thermistors from the same lot were tested, namely: NTC1, NTC2, and NTC3. Figure 14a,b show typical experimental results (red, green, and blue curves) and numerical results (solid and dashed black curves) obtained when NTC1 is immersed with zero immersion length (L i = 0.0 mm) in pure glycerol at rest, using two different constant excitation voltages, v 0 = 7.63 V and v 0 = 6.70 V, respectively. A working temperature T 0 = 25.0 • C was considered here. As shown by the dashed black curves, the ideal model (see Figure 11, where R L = ∞) does not allow for adequate description of the self-heating of the active core (measurements curves 1 to 3), regardless of the excitation voltage used. In contrast, the introduction of a finite loss resistor in the core model (see Figure 12, where R L < ∞) has systematically allowed a very adequate and realistic description of the thermistor active core self-heating, as shown by the black solid curves, which well describe the experimental evolutions of the temperature T c (t).
Thus, it can be seen from the curves of Figure 14 that introducing a finite loss resistance R L in the 1D systemic model of the NTC is crucial to properly describe the evolution of the core temperature as a function of time.  Figure 15a shows the influence of the immersion length L i (and, thus, of the loss resistance R L ) on the self-heating of the thermistor NTC3 when immersed in glycerol at rest, at the working temperature T 0 = 25.0 • C. The values of the loss resistance corresponding to each length L i (see Figure 15b) and the value of the contact resistance R c (see Table 3) were determined by minimization of (21), using the same number of samples N f = 300 for each set of measurements. From the modeling point of view, both the equivalent insulating sheath and the fluid domain were discretized with the numerical scheme exposed in Section 3.1, using the same number M of layers (M = 28) in both cases. The radius a f of the fluid domain was chosen equal to that of the measurement cell: a f = 4.5 mm. The values of the other parameters required for the numerical modeling (except R L ) were those indicated in Table 3. As it can be seen from the curves shown in Figure 15a, the self-heating of the thermistor core is greater the shorter the immersion length, all other factors being equal. This clearly indicates that there are heat exchanges between the connecting wires and the liquid to be characterized and that these exchanges are greater the longer the length of wire immersed, as expected. However, the non-linearity of the variations of R L with L i (see Figure 15b) reflects the great complexity of the thermal exchanges between the connecting wires and the fluid to be characterized. Moreover, the fact that R L is not infinite when L i = 0.0 mm also reflects the existence of thermal transfers between the NTC core and the connector system which connects the electrical control circuit to the NTC copper wires, which are acting as thermal paths.
The heat sink property of the connecting wires is also well illustrated by the results gathered in Table 4. When L i > 0.0 mm, it can be seen from this table that the thermal power current i th L , which is evacuated via the connection wires, represents more than 50% of the total thermal power i tḣ q e supplied by the control circuit to the NTC at t = t f and is, therefore, not used for the thermal characterization of the fluid. However, it can also be observed from this table that the thermal power fraction i th C , that only serves to raise the active core temperature, has become negligible at t = t f , whatever the value of L i . Meanwhile, the fraction i th S+F+T 0 of the thermal power that passes through the sheath, then the fluid, and, finally, reaches the thermostat is considerable and allows, despite the losses, for adequately probing the fluid to be characterized, as well as for the four L i values considered in this study.  Table 4 confirm the relevance of the modeling proposed in this paper.
In conclusion, we found in this study that the value of the loss resistance R L considered in the systemic model was highly dependent on the immersion length L i . This dependence has a considerable influence on the results provided by the modeling and, thus, on the good correspondence between numerical and experimental results (see Figure 14a,b, for example). In the context of the present approach, therefore, it is essential to access to the value of the loss resistance with a high degree of accuracy. Unfortunately, the value of the loss resistance corresponding to a given immersion length L i = 0 depends on the nature of the considered liquid to be characterized. Therefore, when L i = 0, it is impossible to determine both R L and another unknown characteristic of the fluid (such as thermal conductivity, for example). Therefore, it is generally recommended, on the one hand, to limit as much as possible thermal exchanges between the connection wires and the liquid to be characterized (for example, by limiting the immersion length of the connecting wires or by highly isolating them) and, on the other hand, to ensure a very precise and reproducible positioning of the NTC (because of the high value of the slope of R L (L i ) when L i → 0).
Among the four immersion lengths considered in this study, we noted that the length L i = 0.0 mm was the most relevant, which is consistent with the previous discussion. Indeed, it can be seen from Table 4 that this specific length allows, on the one hand, for obtaining the lowest thermal losses; on the other hand, it leads to the highest fluid probing thermal current i th S+F+T 0 , which is in favor of good thermal characterizations. Moreover, the heat exchanges between the fluid to be characterized and the connecting wires are negligible in the case where L i = 0.0 mm. Therefore, the value of the loss resistance becomes fluid independent in this particular case, which is a great advantage for thermal characterization, as mentioned previously.

Thermal Power Balance and Characteristics Times
The 1D systemic modeling introduced in this paper also allows access to quantities that are difficult or impossible to measure, such as the various thermal current densities that flow through the system during thermal characterization.
The analysis of the temporal variations of the thermal powers exchanged during the thermal characterization of glycerol provides useful information about the operation of a NTC in thermal characterization mode. Figure 16 shows the evolution of the thermal power current densities (in W/m 2 ) calculated at position r = a c as a function of time t ∈ [0, 30] s.  Analyzing the curves shown in Figure 16 reveals the existence of mainly three different temporal domains, delimited by two characteristics times t c1 and t c2 for which we get i th , respectively. The operation of the NTC on each of these three temporal domains can be described as follows: • when 0 t < t c1 , time domain (I): only the thermal current density i th C is significant, and the thermal power supplied by the electrical control circuit mainly serves here to raise the temperature T c of the NTC active core. This time domain cannot be used to characterize the fluid surrounding the NTC because the evolution of T c (t) between 0 and t c1 is mainly influenced here by the physical properties of the active core and not by the physical properties of the surrounding fluid to be characterized, that is not yet probed by the thermal waves emitted by the NTC. In addition, we can see in Figure 15a that the curves giving the experimental temporal evolution of T c coincide with t ∈ [0, t c1 ) and this, whatever the value of L i and, thus, of R L . This observation confirms that the evolution of T c (t) over the interval [0, t c1 ) mainly reflects the physical properties of the thermistor active core only. Therefore, it is this particular time domain that was used in the present systemic modeling to determine the optimal values of the effective a c and ρ c c c parameters, by minimizing (21) on [0, t c1 ) (see Table 3). • When t c1 < t < t c2 , time domain (II): in the case of Figure 16, the three thermal current densities i th C , i th L , and i th S+F+T 0 have comparable values; therefore, the thermal characterization of the surrounding fluid using the time evolution of the core temperature will be influenced by both the properties of the core, of the insulating shell and by the thermal losses via the connecting wires. Given these various influences, the usual transient methods of thermal characterization (which do not use systemic modeling) based on bead-type NTCs should avoid the use of domain (II) data. • When t > t c2 , time domain (III): in this case, the thermal current density i th C has become negligible in front of i th L and i th S+F+T 0 , and it can be supposed that the properties of the core are then without any significant influence on the temporal evolution of its temperature T c . Consequently, it is this time domain that should be used preferably for thermal characterization of materials by using bead-type NTCs self-heating methods. However, since thermal losses and the insulating sheath still have an influence on the temporal evolution of T c in the (III) domain, thermal characterization methods that do not use systemic modeling must imperatively resort to a prior calibration of the measuremen device by using several reference fluids, such as glycerol, ethanol, water-glycerol mixtures, and gelled water (using agar-agar, for example). Note that, in this work, the characteristic time t c2 has been calculated using the following relation: i th The 1D systemic modeling proposed in this work allowed to highlight the existence of two characteristic times of the thermal transfers, when bead-type NTCs are used for thermal characterization of materials by self-heating methods. The knowledge of these characteristic times can help to a better understanding of the operation of bead-type NTCs used for thermal characterization of materials, as it will be illustrated in the following paragraph.

Constant Voltage Excitation Signal Processing
Various kinds of possible electrical excitations can be consider to heat the NTCs used in transient thermal characterization methods: current step through the NTC (using a precise current source) and resistance step (used in constant temperature heating techniques (CTHT)), and voltage step (using a divider bridge). From the point of view of the control circuit, thermal characterization using a constant voltage excitation is relatively easier to implement than CTHT; thus, it is very common to use a voltage step to heat NTCs in transient thermal characterization. Several signal processing methods can be considered in order to extract the thermal conductivity from the temporal variations of T c (t) in the case of a voltage step excitation.
For example, Kharalkar et al. proposed to write the quantity P e (t)/δT c (t) in the following form [11]: They then proposed to express the thermal conductivity k f of the fluid to be characterized as k f = 1/(a 1 D −1 0 + a 2 ), where a 1 and a 2 are empirical coefficients depending on both the thermistor considered and the thermal losses. These coefficients are determined by a precise calibration, at given voltage excitation, temperature, and immersion length, using at least two reference materials with thermal conductivities close to the one being measured.
In the present paragraph, we consider the processing of the experimental data shown in Figure 15a, obtained when L i = 5.0 mm. As can be seen from Figure 17a  Other signal processing can also be considered. As can be seen from Figure 17b, a modeling in the form of a CTHT type law P e (t) = P e,β t −1/2 + P e,Γ also provides an acceptable description of the experimental measurements, even in the case of constant voltage heating, provided only data of time domain (III) are considered.

Liquids Thermal Conductivity Measurements
In this paragraph, a new method of thermal conductivity measuring is introduced, which we describe as quasi-absolute, since it requires the use of only one reference liquid.

Electro-Thermal Systemic Modeling (ESM) Method
The present 1D electro-thermal systemic modeling (ESM) approach was applied to the determination of the thermal conductivity k f of three different water/glycerol mixtures, using the following water mass ratios: 100% water (100W0G), 50% water (50W50G), and 40% water (40W60G). Thermistor NTC3 was used in constant voltage heating mode (using the circuit of Figure 4), with L i = 0.0 mm.
The implemented protocol was as follows:  Tables 3 and 4) were used in this work for the determination of the thermal conductivities of the three water/glycerol mixtures.

2.
Voltage step excitation of the NTC (using the circuit of Figure 4) precisely immersed at L i = 0.0 mm in the liquid to characterize, at constant working temperature, T 0 . The experimental time variations of the NTC core temperature T mea c (t) were extracted from the v 0 (t) and v 1 (t) voltages, recorded using a data acquisition board (see Section 2.3).

3.
Determination of the measured thermal conductivity value k mea f by minimization of δT rms c , given by Equation (21), as a function of the thermal conductivity value k f used in the systemic model. The values of the fluid density ρ f and its specific heat c f are supposed known.

Pure Water Liquid (100W0G)
In the case of pure water (100W0G), a small amount of agar-agar powder has been added to 100W0G samples in order to limit the influence of free convection. Figure 18a shows typical experimental (blue line) and 1D ESM computed (black and red lines) time evolutions of the NTC active core temperature δT c = T c − T 0 , when the thermistor is immersed at L i = 0.0 mm in pure water at rest, at a working temperature T 0 = 24.0 • C and using a voltage step excitation v 0 = 8.46 V. The electrical circuit used to both excite the NTC and record the data is that shown in Figure 4, using the same settings as those chosen for the thermal characterization of glycerol (see Section 4.2.1). The correspondence between the experimental measurements (blue curve) and systemic modeling (black curve) is very good when the value of k f is set in the model to k mea f = 0.586 W · K −1 · m −1 , as reflected by the low value of δT rms c in this case. In addition, note that the very good correspondence between calculated and experimental curves over the time domain [0, t c1 ] confirms the relevance of the NTC core effective physical parameters values obtained from the thermal characterization of pure glycerol.   The red dashed curves shown in Figure 18a illustrate the δT rms c minimization steps. By gradually varying the value of k f used in the systemic model (from k f = 0.5 to 0.7 W · K −1 · m −1 , for example), all other quantities being held constant, δT rms c is varied accordingly (see Figure 19). Thus, the value of k f that gives the smallest value of δT rms  Table 1.
Finally, the curves presented in Figure 18b show in the present configuration (pure liquid water at rest, v 0 = 8.46 V, and L i = 0.0 mm), when t > t c2 , that the thermal current density i th S+F+T 0 is quite larger here than both i th L and i th C . Thus, it can be concluded that bead type NTC self-heating methods are in the present configuration quite adequate for thermal characterization of pure water at rest.

Glycerol-Water Mixtures 50W50G and 40W60G
The same measurements as in Section 4.3.2 were reproduced both in the case of 40W60G and 50W50G mixtures, under the same conditions as in the case of pure water. The electrical circuit used to both excite the NTC and record the experimental data was that shown in Figure 4, using the same settings as those chosen for the thermal characterization of pure glycerol at rest.
The values of the glycerol-water mixtures densities ρ f and specific heats c f used in the 1D systemic modeling were taken from the open-source thermophysical properties library CoolProp [40].
Glycerol-Water Mixture 40W60G Figure 20a shows typical experimental (blue line) and ESM computed (black lines) time evolutions of the NTC active core temperature δT c , when the thermistor is immersed at L i = 0.0 mm in 40W60G mixture at rest, at a working temperature T 0 = 24.0 • C and using a voltage step excitation v 0 = 8.39 V. The correspondence between the experimental measurements (blue curve) and the electro-thermal systemic modeling (black curve) is very good here when the value of k f is set to k mea f = 0.384 W · K −1 · m −1 in the model, as reflected by the low value of δT rms c obtained in the present case. Finally, as in the case of pure water thermal characterization, the curves presented in Figure 20b show that the thermal current density i th S+F+T 0 is quite larger here than both i th L and i th C when t > t c2 . Therefore, it can, thus, be concluded that methods based on self-heating of bead-type NTCs are quite adequate for thermal characterization of 40W60G mixtures.   As in the case of pure water and 40W60G mixture thermal characterizations, the curves presented in Figure 21b show that the thermal current density i th S+F+T 0 is again quite larger here than both i th L and i th C . Therefore, it can still be concluded here that methods based on self-heating of bead-type NTCs are quite adequate for thermal characterization of 50W50G mixtures.

Synthesis
The series of measurements exposed in the previous paragraphs (Sections 4.3.2 and 4.3.3) were repeated N = 12 times for each of the three liquids studied in the present work. A statistical analysis was then applied to these measurements, leading to the results collected in Table 5. Table 5. Determination of the thermal conductivities of several liquids (pure water and waterglycerol mixtures) by the present 1D ESM approach, with: k the mean value of measurements, δk = σ k / √ N − 1 the standard error of measurement, σ k the standard deviation of measurements, and k ref the reference thermal conductivity value (all these quantities are given in W · K −1 · m −1 ). The results gathered in Table 5 show that the 1D ESM approach introduced in this paper gives very suitable results for each of the three liquids that were characterized in the present study. Although this new approach uses only one reference liquid (pure glycerol here), it can be seen from Table 5 that this method allows acceptable measurements of thermal conductivities of liquids with large thermal contrasts (e.g., pure water and 40W60G mixture), as well as liquids with fairly close thermal conductivities (e.g., 40W60G and 50W50G mixtures).

Liquid
However, it is also observed from Table 5 that the standard error of measurement δk decreases with increasing glycerol content. This decrease is probably due to the fact that it is glycerol that has been used here as the reference liquid. Despite this tendency, it was found in all cases that the measurements were quite accurate. Even in the worst case (pure water here), where the value δk/k = 2.3 % was obtained, it was found that the measured value k = 0.607 W · K −1 · m −1 was very close to the reference value k ref = 0.605 W · K −1 · m −1 .

Concluding Remarks and Perspectives
A new method for measuring the thermal conductivity of liquids, with the advantage of requiring only one reference liquid, has been introduced and described in this work. This new quasi-absolute approach has proved to be easy to implement, fast (each series of measurements typically takes less than one minute), and quite accurate, regardless of the thermal conductivity range covered during the measurements.
The development of this new thermal characterization method was made possible through the use of a realistic 1D electro-thermal systemic modeling of the whole system (Fluid/Sensor/Connectics/Control electrical circuit), based on a Godunov-SPICE type discretization scheme. This systemic modeling has allowed to highlight two characteristic times of heat transfers within a self-heated thermistor, allowing for better understanding of the operating of bead-type NTCs when they are used for thermal characterization of liquids at rest.
The discretization scheme and the 1D systemic modeling (including the NTC sensor, the electrical control circuit and the materials to characterize) used in the present study were fully described and documented in the paper. Note that the different computer codes (written in Python 3.8) used for the control of the electrical circuit and the processing of the experimental measurements, as well as the SPICE scripts used for computing the systemic modeling numerical results, are available upon request from the author.
The perspectives of this work are numerous. First, from the point of view of thermal characterization, other viscous reference materials (such as pastes and gels) will be tested, and the method will be applied to the measurement of a wide range of thermal conductivities, as a function of temperature. We believe that the precise positioning of the NTC sensor at the specific immersion length L i = 0.0 mm can quite easily be automated. This precise and automated positioning of the sensor could allow the present method to access a wide variety of fields, such as the thermal characterization of biological materials (in-vitro and in-vivo) and nanofluids, in the laboratory or in industrial processes. We believe that this approach is also of interest in the case of dynamic processes, such as the drying of complex suspensions, where the thermal conductivity and the level of the free surface of the suspension can change over time.
By changing the boundary condition considered in the present 1D systemic modeling (at the interface between the insulating shell and the fluid), it will also be possible to consider the application of the present approach to the measurement of fluids thermal characteristics when they are flowing.
Systemic 1D modeling of positive temperature coefficient (PTC) silicon thermistors (silistor) is also considered as a possible extension of the present work.

Data Availability Statement:
The data supporting the reported results were conducted by the author and can be found from him.