*2.1. Experimental Set-Up*

The experimental heating tests consisted of heating a vessel to 200 ◦C for 1800 s on an induction cooktop prototype from BSH Home Appliances. The inductor used in the experiments generates a power distribution resembling a ring, which is assumed to be rotationally symmetric; see Figure 1 in Cabeza-Gil et al. [10]. The applied power turns into heat in a steel microlayer of 100 μm, which is placed at the bottom of the vessel, by means of the dissipation of the eddy current density induced [11]. The power density was controlled with a PI algorithm, whose input is a thermocouple, hereinafter referred to as *Tsensor*, which was located at the centre of the cookware surface for all vessels [12]. To measure the thermal footprint in the glass, four thermocouples were placed under the glass at a radius of 7.00, 5.00, 2.00 cm and at the centre, *T*1, *T*2, *T*<sup>3</sup> and *T*4, respectively, see Figure 1b.

Three different solids (two vessels and a steel plate) with different concavities, see Figure 2, were heated: a multilayer WMF pan with a diameter of 21 cm and a thickness of 4.9 mm (red line); a multilayer Schulte pan, model Industar, with a diameter of 20 cm and a thickness of 3.9 mm (green line); and a circular steel plate, which was specifically used in this study since its base surface is practically flat, with a diameter of 20 cm and a thickness of 6 mm (blue line). The multilayer pans consisted of three layers: steel, aluminium and steel, from bottom to top. The WMF has an aluminium volume percentage of 71.43%, while the Schulte pan has an aluminium volume percentage of 74.36%. The concavity of each vessel was measured with a Faro Prime robot with an accuracy of ±27.00 μm (see Figure 2). The emissivity () of each vessel under investigation was measured with a thermal emissometer, model TIR 100-2 from Inglas, along the whole spectrum.

**Figure 2.** The mean and deviation concavity for each solid under investigation are depicted. The red line corresponds to the WMF, the green line to the Schulte and the blue line to the steel plate. The axis x = 0 is placed at the lowest geometric point of the vessel located in the periphery.

A computer to which the experimental setup was connected calculated the supplied power required and provided it to achieve the target temperature (200 ◦C) by means of the PI controller [12]. The tests for each sample were carried out three times. The software used to supply the power and to register all temperatures of the thermocouples was MATLAB R2020a.

#### *2.2. Model Description*

#### 2.2.1. Heat Transfer Model

Two different FE models were developed to simulate the heating process of the vessels depending on how the heat transfer between the vessel and the glass was modelled. Both models were divided into two coupled systems, the vessel (*S*1) and the glass surface (*S*2). A new solid domain, a layer of air (*S*3), was included in Model I to simulate the thermal contact between the vessel and the glass. In Model II, this thermal contact was modelled through a variable thermal contact conductance.

In *S*1, *Q*˙ *pan* represents the heat rate generated in the ferromagnetic microlayer of the pan from the induction cooking hob. *<sup>Q</sup>*˙ *pan*−*amb* indicates the convective and radiative heat losses from the pan to the ambient environment, modelled by *h*<sup>1</sup> and *hr*1, respectively, which are the heat transfer convective and radiation coefficients. Heat losses between the pan and the glass are different in each model because of the presence of air. In Model I,

*<sup>Q</sup>*˙ *pan*−*air* depicts the heat losses to the air between the pan and the glass. In *<sup>S</sup>*2, *<sup>Q</sup>*˙ *air*−*gla* represents the conductive heat transfer between the air and the glass and *<sup>Q</sup>*˙ *pan*−*gla* between the pan and the glass. In Model II, without the layer of air, *<sup>Q</sup>*˙ *pan*−*gla* depicts the conductive heat transfer between the pan and the glass. *Q*˙ *gla* is the heat absorbed by the glass, and *<sup>Q</sup>*˙ *gla*−*amb* are the convective and radiative losses from the glass to the ambient environment, which are also modelled as *h*<sup>1</sup> and *hr*1, respectively.

The governing equations of the systems are described in Cabeza-Gil et al. [10], which come from the local heat transfer equation (Equation (1)) [13,14].

$$\begin{cases} \quad (a) \ P = \rho\_{SM} c\_{\text{c-SM}} \frac{\partial T\_{\text{SM}}}{\partial t} - k\_{SM} \nabla^2 T\_{SM\prime} \Rightarrow \text{ SM}, \\\\ \quad (b) \ 0 = \rho\_{SD} c\_{\text{c-SD}} \frac{\partial T\_{\text{SD}}}{\partial t} - k\_{SD} \nabla^2 T\_{\text{SD}\prime} \Rightarrow \text{SD}. \end{cases} \tag{1}$$

The domain of Equation (1a) is the steel microlayer (SM) of the vessel where the electromagnetic power is supplied, being *P* the volumetric power density generated by the induction heat source. Equation (1b) refers to the remaining part of the solid domain (SD). The terms *ρSD*, *ce*−*SD*, and *kSD* are the density, the specific heat capacity and the thermal conductivity of the solid material (steel or aluminium), respectively. *TSM* and *TSD* are the corresponding temperature at some determined point in the volume domain.

Regarding the boundary conditions of the system, Equations (2) and (3) refer to the convection and radiation heat losses to the ambient environment, respectively.

$$1 - \lambda\_i \frac{\partial T\_i}{\partial n} = h\_i (T\_i - T\_{amb}) \, . \tag{2}$$

$$h\_i = h\_i^{conv} + h\_i^r = h\_i^{conv} + \sigma \epsilon (T\_i^2 + T\_{amb}^2)(T\_i + T\_{amb}) \tag{3}$$

where the subscript *i* refers to outer surface of the vessel or the glass and *<sup>∂</sup>Ti <sup>∂</sup><sup>n</sup>* is the partial T-derivative normal to the reference surface. *hi* includes both the convective and radiative contributions, i.e., *hconv <sup>i</sup>* and *<sup>h</sup><sup>r</sup> <sup>i</sup>* are the convective and radiative heat transfer coefficients, respectively. *σ* is the Stefan–Boltzmann constant, is the emissivity of the pan and *Tamb* is the ambient temperature.

Thermal conduction was differently modelled for Model I and Model II. In Model I, the conduction heat losses were modelled through Equation (4) from the vessel to the air (*hva <sup>c</sup>* ) and from the air to the glass (*h ag <sup>c</sup>* ):

$$-\lambda\_v \frac{\partial T\_v}{\partial \mathbf{n}} = h\_c^{va} (T\_v - T\_a)\_\prime \ -\lambda\_a \frac{\partial T\_a}{\partial \mathbf{n}} = h\_c^{\mathbf{g}\mathbf{g}} \left(T\_a - T\_{\mathbf{g}}\right) \tag{4}$$

where *λ<sup>v</sup>* and *λ<sup>a</sup>* are the thermal conductivity of vessel and air. *Tv*, *Ta* and *Tg* are the respective temperatures at the vessel, air and glass surfaces. *hva <sup>c</sup>* is the thermal contact conductance to be evaluated between the vessel and the air, and *h ag <sup>c</sup>* is the thermal conductance between the air and the glass. A perfect thermal contact between the vessel and the air surface, and the air and the glass surface, was considered.

For the Model II, an iterative analysis was performed in the three samples under investigation to determine the relationship between the air gap and the thermal conductance coefficient along the radius, see Equation (5):

$$Q\_{\mathcal{L}} = \int\_0^R \int\_0^{2\pi} h\_{\mathcal{L}}^{\mathbb{R}\mathbb{Z}}(r) \cdot \Delta T(r) \cdot r \cdot dr \cdot d\theta \,, \tag{5}$$

where *Qc* refers to the heat loss at the contact between the pan and the glass, *h vg <sup>c</sup>* (*r*) corresponds to the thermal conductance coefficient depending on the radius between the vessel and the glass, and Δ*T*(*r*) is the difference in the surface temperatures in the volume domain.

#### 2.2.2. Finite Element Model

The FE model developed in this study is shown in Figure 3, and it is based on the study developed by Cabeza-Gil et al. [10]. The vessel geometry and the inclusion of the layer of air were different depending on the FE model analysed. Model I includes an air layer between the solid and the glass whereas Model II includes a variable thermal contact conductance. Both FE models consisted of a vitroceramic circular glass with a 200 mm radius × 4 mm thickness and the corresponding vessel previously described in Section 2.1. Due to the rotational symmetry, a quarter of the model was designed.

**Figure 3.** Vitroceramic glass is represented in orange, and the multilayer round WMF pan is represented in blue. Model I is modelled with the air (purple instance) between the glass and the pan, while Model II replaces the air by modelling the interaction between the vessel and glass by a variable thermal conductance coefficient along the radius.

The commercial software Abaqus v.6.14 was used to generate the model mesh and perform the simulations. The glass was considered a shell, which was meshed with 1408 quadratic 4-node quadrilateral (DS4) and 32 3-node linear triangular (DS3) heat transfer shell elements, while the vessel (WMF, Schulte pan and steel plate) was considered a 3D solid, which was meshed with approximately 35,000 (depending on the solid) 10-node quadratic heat transfer tetrahedral elements (DC3D10); see Figure 3. For the Model I, the layer of air was approximately modelled with 200 DC3D10 elements.

The power density distribution generated by the induction system was numerically computed in an electromagnetic FE model [15,16]. The power-density field from the FE model was mapped using an in-house subroutine, which was written in MATLAB 2020a, onto the FE mesh to perform the thermal analysis. The power supplied was calculated using PI control, which is a temperature-level control, and it was added through the URDFIL and DFLUX subroutines. This control reproduces the control implemented in the experimental heating tests by determining the power supplied from the previous time increment (Δ*t* = 1 s) of the temperature sensor. The target temperature in the control was 200 ◦C, and the maximum power applied was limited to 2200 W, as in the experimental tests.

The whole model was at ambient temperature (*Tamb* = 23 ◦C) as initial condition. The boundaries conditions applied to the model are described in Section 2.2.1. Briefly, convection and radiation heat losses were imposed to all external surfaces of the model. For Model I, the heat conduction between the vessel and the air, and the air and the glass, was considered through a thermal contact conductance. Whereas for Model II, the heat conduction between the vessel and the glass surfaces was considered through a variable thermal contact conductance depending on the position (*h vg <sup>c</sup>* (*r*). Convective and conductance parameters used in the computational models were obtained individually for each vessel after an optimisation process by fitting the experimental tests with the numerical tests. In each vessel, the same convective parameter was applied for all external surfaces. The coefficients of the glass were the same for all simulations. The thermal air properties [14] used in the study are shown in Table 1.


**Table 1.** Properties of air (Model I) modelled based on temperature as a continuum between the vessel and the pan [12].

2.2.3. Finite Element Method Discretization

The weak form of Equation (1a) can be written as:

$$\int\_{\Omega} P\delta T d\Omega = \int\_{\Omega} \rho c\_{\varepsilon} \frac{\partial T}{\partial t} \delta T d\Omega - \int\_{\Omega} k \nabla^2 T \delta T d\Omega \tag{6}$$

being Ω the solid domain and *δT* the virtual temperature. The application of classical differentiation rules to the last term of Equation (6) leads to the following statement:

$$\int\_{\Omega} P\delta T d\Omega = \int\_{\Omega} \rho c\_{\epsilon} \frac{\partial T}{\partial t} \delta T \, d\Omega + \int\_{\Omega} \nabla \delta T \, k \nabla T d\Omega - \int\_{\partial \Omega} k \nabla T \delta T \, d(\partial \Omega) \tag{7}$$

The FEM discretization procedure starts from the following approximation of the temperature function:

$$T(x, y, z) = N\_a T\_{a\prime} \text{ being } \alpha = 1, \dots, n \tag{8}$$

where *T*(*x*, *y*, *z*), the temperature in the Cartesian coordinates, is represented by *Nα*, the shape functions, and *Tα*, the nodal temperatures. *n* is the total number of degrees of freedom of the model. Following the usual approximation of FEM, the virtual temperature (*δT*) is identified as the shape functions, *δT* = *Nβ*. The matrix form of Equation (7) is represented as follows:

$$Q\_a^{(\varepsilon)} = \int\_{\Omega^{(\varepsilon)}} P \mathcal{N}\_{\hat{\beta}} d\Omega + \int\_{\partial \Omega^{(\varepsilon)}} k \nabla T \mathcal{N}\_{\hat{\beta}} \, d(\partial \Omega) \tag{9}$$

$$K\_{a\mathcal{S}}^{(\varepsilon)} = \int\_{\Omega^{(\varepsilon)}} (k\_x N\_{a,x} N\_{\mathfrak{F},x} + k\_y N\_{a,y} N\_{\mathfrak{F},y} + k\_z N\_{a,z} N\_{\mathfrak{F},z}) d\Omega \tag{10}$$

$$M\_{\mathbf{a}\boldsymbol{\beta}}^{(\varepsilon)} = \int\_{\Omega^{(\varepsilon)}} \rho c\_{\varepsilon} N\_{\mathbf{a}} N\_{\boldsymbol{\beta}} d\Omega \tag{11}$$

The assemblage and condensation procedures for the system matrices and vectors lead to the well-known final system of algebraic equations, which has been solved implicitly through the trapezoidal rule for time integration in Abaqus. The initial and boundary conditions needed to solve the equations system are described in Section 2.2.1.

$$M\_{a\beta}\frac{\partial T\_{\beta}}{\partial t} + K\_{a\beta}T\_{\beta} = Q\_a \tag{12}$$

#### *2.3. Design of the Experiment*

The influence of the key parameters of the model, the conductivity (*k*), specific heat (*ce*), emissivity (), and concavity (*Con*) of the vessel, and the convective coefficients of the vessel (*hconv <sup>v</sup>* ) and the glass (*hconv <sup>g</sup>* ) in the cooking process were analysed following the DoE methodology by a full factorial analysis [17]. The simulations consisted in heating the vessel during 1800 s as in the experimental tests. A screening analysis was performed to observe the influence of each variable and decide the levels of the DoE. Based on this analysis, an intermediate value was selected for the conductivity, specific heat and concavity, whereas the remaining terms (emissivity and convective coefficients of the vessel and the glass) had two levels, resulting in 243 simulations (i.e., 33 × 32 = 243 simulations), see Table 2.

**Table 2.** Values of the analysed parameters for each level. Conductivity, specific heat and concavity have three levels, whereas both convective parameters and emissivity have only two. The steel density was considered constant with a value of 7900 kg/m3.


The maximum and minimum levels of the conductivity and the specific heat were chosen based on the properties of aluminium and steel [18]. The parameters affecting three heat losses (, *hv*, *hg*) were chosen according to values reported in the literature [10,11].

To analyse the influence of concavity, three different scenarios, see Figure 4, were considered. The red line was considered the reference case (from the WMF pan), hereinafter referred to as scenario #B. Two more scenarios where the concavity was increased by 1.5 times and decreased by 0.5 times were introduced (scenario #A and scenario #C, respectively).

**Figure 4.** Concavities of the three different scenarios. Scenario #A (blue line) corresponds to a larger concavity compared to scenario #B, the concavity measured to the WMF pan with the robot (red line). Scenario #C (green line) corresponds to a flatter concavity compared to scenario #B.

The effect of the key parameters in the heating process were analysed measuring different factors of the cooking process: time to reach steady state (*tst*), supplied energy during the cooking (*Ein*), maximum temperature of the sensor (*Ts*) and temperature homogenisation at the bottom along the radius after 1 min from the start (*Hr*):

$$H\_r = 1 - \frac{1}{S\_p} \int\_0^{2\pi} \int\_0^{r\_{ext}} \frac{|\,\,\tilde{T} - T(r, \theta)\,\,|}{\tilde{T}} \cdot r \cdot dr \cdot d\theta \,\,\tag{13}$$

where *Hr* refers to the temperature homogenisation along the radius in the second 60 of the cooking, *Sp* is the cookware surface of the pan, *r* is the radius of the pan, *T*¯ is the mean temperature of the nodes selected in the cookware surface and *T*(*r*, *θ*) is the temperature in each node.

#### **3. Results & Discussion**

### *3.1. Experimental-Numerical Calibration of the Heat Loss Coefficients*

The heat loss coefficients for the five calibrated models are presented in Table 3. Only the WMF and Schulte pans were developed in Model I because the steel plate is flat (Figure 2) and there is no layer of air to be modelled. Convective coefficients were optimised to reduce the mean absolute error (MAE) between the experimental and numerical temperatures. We obtained a convective coefficient of 8.0 W/m2K, similar to Sanz-Serrano et al. [11], which obtained a convective coefficient of 9.5 W/m2K using the difference finite method.

**Table 3.** Heat loss coefficients: convective coefficients (*hconv <sup>v</sup>* and *hconv <sup>g</sup>* ), thermal conductance coefficients (*h vg <sup>c</sup>* , *hva <sup>c</sup>* and *h ag <sup>c</sup>* ) and emissivity coefficient () for the five FE models, the two with the layer of air (Model I) and the three with the non-linear thermal conductance along the radius (Model II).


3.1.1. Model I: Modelling the Layer of Air between the Vessel and the Pan

Figure 5a,b show the WMF and Shulte heating processes, respectively. The dotted lines indicate the values obtained from the experimental tests whilst the continuous lines represent the computational results (see Figure 1b for the location of the thermal sensors). Shaded areas in the experimental tests are the standard deviations. Both numerical results of the surface thermocouples are within the experimental deviation.

**Figure 5.** Temperatures of the surface and area under the glass from 0 to 1800 s for both the experimental test (dotted lines) and computational simulations (continuous lines) in Model I. (**a**) corresponds to the WMF pan and (**b**) to the Schulte pan.

The fitting between experimental and numerical results is good, although the experimental glass temperatures of the WMF pan heat up slightly faster than the simulated temperatures in the transient state (0–200 s). The temperatures in the steady state are

practically the same. Table 4 shows the MAEs between the numerical and experimental temperatures. The thermocouple placed in the centre (*T*4) presented the highest error, 6.75 ◦C.

**Table 4.** MAE of the three models with Model I. *Tsensor*, *T*1, *T*2, *T*<sup>3</sup> and *T*<sup>4</sup> are the thermocouples in Figure 1b. The units are ◦C.


On the other hand, during the whole heating process in the Schulte (Figure 5b), both the experimental and simulated results are similar (all MAEs are under 5 ◦C) except for *T*3, which differs more from the experimental temperature (MAE of 13.27 ◦C). A misalignment of the thermal sensor or the vessel in the experimental tests might be one factor. This area has the largest concavity deviation and small variations can significantly influence the thermal contact conductance and thus the temperature.

#### 3.1.2. Model II: Fitting of the Nonlinear Thermal Conductance along the Radius

The equation that relates the thermal conductance coefficient and the air gap for the pans (WMF and Schulte) is shown in Figure 6. An iterative inverse analysis, such as in Paesa et al. [19], was performed, and it included several simulations of WMF with different variable thermal conductances along the radius. The optimal result was correlated with the concavity of the pan (see Figure 2), and thus, the following relationship, see Figure 6, was achieved. The WMF pan was used as reference vessel and the parameter fitting were also employed later for the Schulte pan simulations.

**Figure 6.** Relation between the thermal conductance coefficient and the air gap. The relation and its equation, where *y* is the distance between the glass and the base of the pan, are shown graphically. The R-squared value of the regression model is also presented.

The results of the simulations are shown in Figure 7. The surface thermocouples of WMF and Schulte present very low MAEs of 3.09 ◦C and 1.80 ◦C, respectively (see Table 5).

**Table 5.** MAE of the two models modelled with Model II. *Tsensor*, *T*1, *T*2, *T*<sup>3</sup> and *T*<sup>4</sup> are the thermocouples in Figure 1b. The units are ◦C.


Regarding the glass thermocouples in the WMF (see Figure 7a), the experimental temperatures are slightly higher than the numerical temperatures in the transient state (approximately 0–200 s). For the Schulte pan, the experimental values are moderately higher than the simulated one during the whole cooking (see Figure 7b).

The heat conduction between the steel plate and the glass was initially considered as perfect as the plate is completely flat; however, the numerical results did not fit the experimental results. Thus, a thermomechanical analysis was performed to observe if there was a significant thermal deformation during the cooking that produced an input of air between the plate and the glass. It was noted that the concavity changes with increasing distance between the base vessel and the glass up to an axial displacement of 120 μm (see Supplemental Data). Thus, *h vg <sup>c</sup>* was modelled as piecewise as follows due to the thermomechanical analysis: from the centre until a radius of 10 cm is 50 W/m2K; from 20 to 50 cm of radii, 150 W/m2K; and the rest of the pan, 500 W/m2K. With this assumption, There was a clear correlation between the experimental and numerical temperatures, see Figure 7c.

**Figure 7.** Temperatures of the surface and under the glass (see Figure 1b for location of the temperature sensors) from 0 to 1800 s for both experimental tests (dotted lines) and computational simulations (continuous lines) in Model II: (**a**) WMF pan, (**b**) Schulte pan and (**c**) steel plate.
