Mathematical Model to Calculate Heat Transfer in Cylindrical Vessels with Temperature-Dependent Materials

Abstract

In this article, a mathematical model capable of simulating the heat transfer of cylindrical vessels whose properties are dependent on temperature is proposed. As a case study, it compares, from an approach of their heat transfer and chemical migration characteristics as a function of the temperature reached, different materials commonly used for the manufacture of water bottles. More specifically, the materials studied were aluminium, polyethylene terephthalate, and polypropylene. The validation of the model consists of an experiment carried out in the laboratory with three water bottles of each of the materials under study, as well as simulations using the Network Simulation Method to recreate the heat transfer that occurs through the walls of the bottles. On the other hand, the nondimensionalization technique is also applied, which allows us to obtain the weight of each of the variables on the problem, as well as the existing relationship between them. Finally, an outside temperature of 30 °C to 50 °C is simulated, which is a common temperature range in southern Europe during the summer season, and an initial temperature of 20 °C for the water contained in the bottle to know the behaviour of the materials and what the final temperature of the water would be after one hour.

1. Introduction

Traditionally in heat transfer problems, the properties of the materials have been considered constant depending on the material to be studied since, in many cases, there was no information available on how they varied as a function of temperature [1,2], but this does not adjust to the real behaviour of the materials, since if we vary the temperature to which the material is subjected, we can observe how the density, thermal conductivity and specific heat, among other properties, its value is altered [3,4,5,6,7,8].

This is why heat transfer problems that are usually carried out considering constant parameters include in their results an error due to the non-consideration of the variability of the properties with the increase or decrease of the temperature. These changes in properties affect the results derived from the heat transfer mechanisms, conduction, radiation and convection [9].

In this work, where the thermal behaviour of a vessel containing a fluid will be studied, as will be shown later, the main mechanisms of heat transfer that can occur are conduction and convection both in its own structure and with its environment, since heat transfer by radiation has been minimised.

The convection process requires the existence of a fluid capable of transporting thermal energy between two points where there is a temperature gradient. Considering that the fluid inside the vessel is at rest, the fluid in charge of transporting the thermal energy will be the moving air surrounding the vessel.

In the radiation process, the existence of an object that transmits heat or direct contact with the heat source is necessary, not generating thermal energy until the electromagnetic (infrared) waves are absorbed. Therefore, in case of application of this study (Section 6), the sun radiates these waves, which are received by the surface of the container which, by absorbing them, transforms them into thermal energy that will subsequently be transmitted to the internal fluid contained in the container through the processes of conduction and convection. For the validation of the proposed mathematical model, radiation to solar exposure will be minimised.

For heat transfer by conduction to exist, there must be a temperature gradient in the volume of a solid or between two materials in contact. In this case, there is a temperature gradient between the surface of the container and the fluid inside it. The transfer of this thermal energy will occur without the exchange of matter due to the vibrational motion of the atoms. As the temperature increases, the atoms will increase their velocity, colliding with neighbouring atoms and transferring this thermal energy to them.

Regarding the properties that influence the conduction mechanism, it is difficult to find conductivity values for many materials, such as some fluids or polymers, and moreover, in many of these cases, they usually present constant values [10]. However, more and more works present this property with temperature dependence [11,12,13,14,15,16], so it is necessary to introduce this variability of the properties with temperature in the mathematical models to improve the accuracy of the simulations.

The proposed mathematical model that will include the aforementioned heat transfer phenomena, as well as the variability of material properties with temperature, is simulated using the Network Simulation Method [9,17,18]. This method has been widely used in numerous problems applied to engineering [19,20], such as diffusion of chlorides in reinforced concrete [21], corrosion problems [22], soil consolidation [23], and solidification problems [24], among others. This model is validated through different experimental tests with three materials, aluminium 319 (Al319), polyethylene terephthalate (PET) and polypropylene (PP). Other authors have used other techniques to solve problems similar to the one posed in this article [25,26].

On the other hand, through the nondimensionalization technique, the relationship between the variables of the problem has been established, allowing us to know the influence of each one of them, as verified in the case study.

As an example of the study of the proposed mathematical model, the distribution of temperatures in bottles of different materials has been studied since the variation of temperatures in plastic materials could lead to an increase in the migration of substances into water [27,28] such as BPA, dimethyl phthalate (DMP), DEHP, OP, among others. Also, this article shows how properties vary as a function of temperature for aluminium 319, polyethylene (PET), and polypropylene (PP), materials commonly used in water bottles. In certain cases, some properties are not significantly affected by temperature, and it can be assumed that not considering the variability would not lead to substantial error, while other properties are affected to a greater extent causing a larger error if considered constant.

In summary, the main objective of this paper is to develop a mathematical model that includes the dependence, in many cases nonlinear, of the material properties on temperature, thus improving the accuracy of the simulations of the network model developed from the mathematical equations. It should be noted the sensitivity of the model to the variability of the properties with temperature is due to the strong coupling of the equations and the nonlinearity, so it is necessary to implement the properties of the materials with the best possible fit. Finally, a study has been added that indicates the weight of each of the variables on the proposed problem carried out using the nondimensionalization technique.

2. Mathematical Model

The mathematical model that is going to be described below is that of a cylinder that can be formed from different materials and whose properties vary with the temperature. As indicated above, the mechanisms to consider are conduction, convection, and radiation. In conduction processes, Fourier’s law is used to determine the heat transfer:

$$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{y}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}=\frac{1}{\mathsf{\alpha }}\frac{\partial \mathrm{T}}{\partial \mathrm{t}}; \mathsf{\alpha }=\frac{\mathrm{k}}{\mathsf{\rho } {\mathrm{c}}_{\mathrm{e}}}$$

(1)

where T is a temperature (K), α is thermal diffusivity (m²/s), k is thermal conductivity (W/mK), c_e is specific heat capacity (J/kg K), ρ is density (kg/m³), t the time (s), and, finally x, y, z are the spatial coordinates. As previously indicated, the thermal conductivity, density, and specific heat capacity, and therefore the thermal diffusivity, will be expressed as a function of temperature.

Because the problem that we are going to develop is a cylindrical vessel, we are interested in defining the problem in cylindrical coordinates. The equation must be defined by the radius, r, the angle, φ, and the height, z, as follows, Figure 1:

$$\frac{1}{\mathrm{r}}\frac{\partial }{\partial \mathrm{r}}\left(\mathrm{k} \mathrm{r}\frac{\partial \mathrm{T}}{\partial \mathrm{r}}\right)+\frac{1}{{\mathrm{r}}^2}\frac{\partial }{\partial \mathsf{\phi }}\left(\mathrm{k}\frac{\partial \mathrm{T}}{\partial \mathsf{\phi }}\right)+\frac{\partial }{\partial \mathrm{z}}\left(\mathrm{k}\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)=\mathsf{\rho } {\mathrm{c}}_{\mathrm{e}}\frac{\partial \mathrm{T}}{\partial \mathrm{t}}$$

(2)

Given the symmetry of the problem, it can be developed in two dimensions, radial and in height, showing the equation to which the network model will be developed, Equation (2):

$$\frac{1}{\mathrm{r}}\frac{\partial }{\partial \mathrm{r}}\left(\mathrm{k} \mathrm{r}\frac{\partial \mathrm{T}}{\partial \mathrm{r}}\right)+\frac{\partial }{\partial \mathrm{z}}\left(\mathrm{k}\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right)=\mathsf{\rho } {\mathrm{c}}_{\mathrm{e}}\frac{\partial \mathrm{T}}{\partial \mathrm{t}}$$

(3)

If we develop the previous equation, it would be given by:

$$\frac{\mathbf{k}}{\mathbf{r}}\frac{\mathbf{\partial }\mathbf{T}}{\mathbf{\partial }\mathbf{r}}+\mathbf{k}\frac{{\mathbf{\partial }}^2\mathbf{T}}{\mathbf{\partial }{\mathbf{r}}^2}+\mathbf{k}\frac{{\mathbf{\partial }}^2\mathbf{T}}{\mathbf{\partial }{\mathbf{z}}^2}=\mathsf{\rho } {\mathbf{c}}_{\mathbf{e}}\frac{\mathbf{\partial }\mathbf{T}}{\mathbf{\partial }\mathbf{t}}$$

(4)

Finally, if we group the terms to facilitate the subsequent development of the network model, the equation would be given as follows:

$$\frac{1}{\mathbf{r}}\frac{\mathbf{k}}{\mathsf{\rho } {\mathbf{c}}_{\mathbf{e}}}\frac{\mathbf{\partial }\mathbf{T}}{\mathbf{\partial }\mathbf{r}}+\frac{\mathbf{k}}{\mathsf{\rho } {\mathbf{c}}_{\mathbf{e}}}\frac{{\mathbf{\partial }}^2\mathbf{T}}{\mathbf{\partial }{\mathbf{r}}^2}+\frac{\mathbf{k}}{\mathsf{\rho } {\mathbf{c}}_{\mathbf{e}}}\frac{{\mathbf{\partial }}^2\mathbf{T}}{\mathbf{\partial }{\mathbf{z}}^2}=\frac{\mathbf{\partial }\mathbf{T}}{\mathbf{\partial }\mathbf{t}}$$

(5)

It should be noted that the grouping of the properties will facilitate its subsequent development in the network model because, as previously indicated, its value depends on the temperature itself.

On the other hand, the implementation of the convection phenomenon is given by Newton’s law of cooling:

$$\frac{\mathrm{d}\mathrm{Q}}{\mathrm{d}\mathrm{t}}=\mathrm{A}\mathrm{h}\left({\mathrm{T}}_{\mathrm{s}}-{\mathrm{T}}_{\mathrm{e}}\right)$$

(6)

where Q is the heat (J), h is the heat transfer coefficient (W/m²K), A is the heat transfer surface area (m²), T_s is the temperature of the solid surface (K), and finally, T_e is the environmental temperature (K).

Finally, the phenomenon of radiation is given by Stefan–Boltzmann’s Law

$$\frac{\mathrm{d}\mathrm{Q}}{\mathrm{d}\mathrm{t}}=\mathsf{\epsilon }\mathsf{\sigma }\left({{\mathrm{T}}_{\mathrm{s}}}^4-{{\mathrm{T}}_{\mathrm{e}}}^4\right)·\mathrm{d}\mathrm{S}$$

(7)

where ε and σ are the emissivity and the Stefan–Boltzmann constant (W/m²K⁴), respectively.

3. Network Model

The Network Simulation Method is a method used to solve various physical-chemical processes that can be defined by a mathematical model. It consists of converting each one of the summands of the equation into an element of an electrical circuit, thus simulating, in this case, a problem of heat transmission through an electrical circuit. That is, an electrical analogy is established between the variables of the problem and the voltage in a central node. This equivalence has been used in heat transfer problems [29,30]; however, the resolution methodology is different. In this case, the proposed method has been widely used satisfactorily in numerous engineering problems [9,17,18].

Therefore, we must proceed to convert Equation (5) into terms that can be used for the Network Simulation Method. The term that appears by multiplying all the summands on the left-hand side of the equation can be expressed as thermal diffusivity, α, which is temperature dependent because its variables are temperature dependent. On the other hand, the terms of the second derivative, after some mathematical operations, will be implemented as two resistors in series.

$$\frac{1}{\mathrm{r}}\mathsf{\alpha }\frac{\partial \mathrm{T}}{\partial \mathrm{r}}+\mathsf{\alpha }\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}+\mathsf{\alpha }\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}=\frac{\partial \mathrm{T}}{\partial \mathrm{t}}$$

(8)

As indicated, in order to implement the second derivatives as resistances, it is necessary that they be free of variable terms. To do this, they are separated into two addends (Equation (9)), where one, its electrical analogy, would be given by a voltage-controlled current generator, and the other, by two resistors in series.

$$\mathsf{\alpha }\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}=\left(\mathsf{\alpha }-1\right)\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}$$

(9)

Performing the same mathematical operation for the term $\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}$, the equation is as follows:

$$\frac{1}{\mathrm{r}}\mathsf{\alpha }\frac{\partial \mathrm{T}}{\partial \mathrm{r}}+\left(\mathsf{\alpha }-1\right)\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{r}}^2}+\left(\mathsf{\alpha }-1\right)\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}+\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}=\frac{\partial \mathrm{T}}{\partial \mathrm{t}}$$

(10)

On the one hand, any second partial derivative free of variable terms is identified in the circuit as two resistors in series. On the other hand, any second partial derivative or first partial derivative with multiplying dependent terms is identified in the circuit as voltage-controlled current sources, and finally, the time derivative represents in the circuit a capacitor. Therefore, for our circuit it would be as follows:

$\frac{1}{\mathrm{r}} \alpha \frac{\partial T}{\partial r},(\alpha -1)\frac{{\partial }^{2}T}{\partial r^2}$, $(\alpha -1)\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}$ are considered voltage-controlled current sources,

$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}$,$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}$ are considered as two resistors in a series,

$\frac{\partial \mathrm{T}}{\partial \mathrm{t}}$ represents a capacitor.

Taking into account that the terms that include α are temperature-dependent functions, voltage-controlled sources will be considered. Prior to the development of the network model, the expressions that include differential equations are spatially discretised by dividing the space into n volume elements according to Figure 2. Thus, the above equations can be expressed in finite differences. Finally, the centres of the cells are balanced with the electrical currents of each of the elements that make up the circuit.

Using the nomenclature of Figure 2, the spatial discretisation of the first and second derivatives is performed by balances between the central nodes and cell edges [30,31]. Thus, the first derivative in finite differences is given by:

$$\frac{\partial \mathrm{T}}{\partial \mathrm{t}}=\frac{{\mathrm{T}}_{\mathrm{i}+1,\mathrm{j}}-{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}}{\mathsf{\Delta }\mathrm{r}}$$

(11)

where $\mathsf{\Delta }\mathrm{r}={\mathrm{r}}_{\mathrm{i}+1,\mathrm{j}}-{\mathrm{r}}_{\mathrm{i}-1,\mathrm{j}}$ and the second derivative by:

$$\frac{{\partial }^2\mathrm{T}}{\partial {\mathrm{z}}^2}=\frac{\frac{{\mathrm{T}}_{\mathrm{i}+1,\mathrm{j}}-{\mathrm{T}}_{\mathrm{i},\mathrm{j}}}{\frac{\mathsf{\Delta }\mathrm{r}}{2}}-\frac{{\mathrm{T}}_{\mathrm{i},\mathrm{j}}-{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}}{\frac{\mathsf{\Delta }\mathrm{r}}{2}}}{\mathsf{\Delta }\mathrm{r}}=\frac{{\mathrm{T}}_{\mathrm{i}+1,\mathrm{j}}+{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}-2{\mathrm{T}}_{\mathrm{i}-1,\mathrm{j}}}{\frac{{\mathsf{\Delta }}^2\mathrm{r}}{2}}$$

(12)

where the expression (12) is the value of the resistance when implementing the electrical device

$$R=\frac{1}{\frac{{\mathsf{\Delta }}^2\mathrm{r}}{2}}$$

(13)

Therefore, the circuit that would represent our problem, considering all of the above, would be the following, Figure 3:

Finally, the boundary conditions, radiation and convection, are implemented by means of a voltage-controlled current generator. The network model is simulated using circuit simulation software such as NgSpice [32].

The NgSpice software (version ngspice-39) [32] uses computational algorithms that allow solving electrical circuits derived from coupled and nonlinear mathematical models such as the one proposed in this article. These algorithms, derived from Nagel’s work [33], include Gear’s time fixation methods [34], trapezoidal integration [35], and the Runge–Kutta algorithm. The efficiency and precision are obtained by reducing the local truncation error and facilitating stability in the convergence of the numerical solution.

4. Nondimensionalization Technique

In this section, a study will be carried out on the behaviour of variables through their appraisal in dimensionless monomials. For this, the nondimensionalization technique will be used, which is a mathematical methodology that allows us to identify the monomials that control the problem posed [36,37,38]. The greatest difficulty that it presents is the correct choice of the reference variables, so it is necessary to have a deep knowledge of the problem posed. In addition, the technique of discrimination must be introduced, which introduces the variability of the properties in the different dimensions of space. Although the properties in the posed problem vary with temperature, it is going to be assumed that they are constant, a decision that we will justify at the end of the section, for the application of the nondimensionalization technique since it facilitates its application and, this study, is only intended to set how variables are grouped.

The first step is to establish the dimensionless variables that, thanks to the references used, will be in the range of [0–1], Equation (14).

$${\mathrm{T}}^{\prime }=\frac{\mathrm{T}-{\mathrm{T}}_{\mathrm{i}}}{{\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{i}}} \to {\mathrm{r}}^{\prime }=\frac{\mathrm{r}}{\mathrm{R}} \to {\mathrm{z}}^{\prime }=\frac{\mathrm{z}}{\mathrm{H}} \to {\mathrm{t}}^{\prime }=\frac{\mathrm{t}}{\mathsf{\tau }}$$

(14)

where T_i is the initial temperature of the vessel, T_e is room temperature, H is the height, and τ is the time at which the problem reaches a steady state. It should be noted that in the case of the radius, as shown in Figure 1, only the vessel structure without the fluid will be studied, R = r₁ − r₂, since it is where the heating of the container with the air is taking place, although the solution obtained will also be valid in the event that the entire cylinder is solid, R = r₁. On the other hand, it should also be indicated that the study is being carried out with symmetry, and therefore, the results obtained will also be valid for the complete cylinder. Finally, the cylinder is considered to be heating since T_e > T_i. In case it is cooling, T_i > T_e, the procedure to be followed would be the same.

Next, the dimensionless variables are introduced in Equation (8) to obtain the dimensionless equation:

$$\left[\frac{\mathsf{\alpha }\left({\mathbf{T}}_{\mathbf{e}}-{\mathbf{T}}_{\mathbf{i}}\right)}{{\mathbf{R}}^2}\right]\frac{1}{{\mathbf{r}}^{\mathbf{\prime }}}\frac{\mathbf{\partial }{\mathbf{T}}^{\mathbf{\prime }}}{\mathbf{\partial }{\mathbf{r}}^{\mathbf{\prime }}}+\left[\frac{\mathsf{\alpha }\left({\mathbf{T}}_{\mathbf{e}}-{\mathbf{T}}_{\mathbf{i}}\right)}{{\mathbf{R}}^2}\right]\frac{{\mathbf{\partial }}^2{\mathbf{T}}^{\mathbf{\prime }}}{\mathbf{\partial }{{\mathbf{r}}^{\mathbf{\prime }}}^2}+\left[\frac{\mathsf{\alpha }\left({\mathbf{T}}_{\mathbf{e}}-{\mathbf{T}}_{\mathbf{i}}\right)}{{\mathbf{H}}^2}\right]\frac{{\mathbf{\partial }}^2{\mathbf{T}}^{\mathbf{\prime }}}{\mathbf{\partial }{{\mathbf{z}}^{\mathbf{\prime }}}^2}=\left[\frac{\left({\mathbf{T}}_{\mathbf{e}}-{\mathbf{T}}_{\mathbf{i}}\right)}{\mathsf{\tau }}\right]\frac{\mathbf{\partial }{\mathbf{T}}^{\mathbf{\prime }}}{\mathbf{\partial }{\mathbf{t}}^{\mathbf{\prime }}}$$

(15)

As expected, the range of variation does not depend on T’ [21]. Assuming that the changes of the derived factors in Equation (15) are of unit order of magnitude, the equation gives three coefficients, indicated in brackets, and these, two dimensionless monomials that must also be of unit order of magnitude, Equation (17).

$$\frac{\mathsf{\alpha }}{{\mathrm{R}}^2} \frac{\mathsf{\alpha }}{{\mathrm{H}}^2} \frac{1}{\mathsf{\tau }}$$

(16)

$${\mathsf{\pi }}_1=\frac{\mathsf{\alpha }\mathsf{\tau }}{{\mathrm{R}}^2} {\mathsf{\pi }}_2=\frac{{\mathrm{R}}^2}{{\mathrm{H}}^2}$$

(17)

where the monomial π₁ represents the relationship of the temperature change with diffusion phenomena and the monomial π₂ the geometric relationship of the dimensions of the cylinder.

Given that the test is realised under conditions where convection phenomena predominate, as previously indicated, if we apply the nondimensionalization technique to Equation (6), and knowing that $\mathrm{d}\mathrm{Q}=\mathsf{\rho }{\mathrm{c}}_{\mathrm{e}}\mathrm{A}\mathrm{R}\mathrm{d}\mathrm{T}$ and $\mathsf{\alpha }=\frac{\mathrm{k}}{\mathsf{\rho }{\mathrm{c}}_{\mathrm{e}}}$ two coefficients are obtained that provide the monomial π₃.

$$\frac{\mathsf{\alpha }\mathrm{h}}{\mathrm{R}\mathrm{k}} \frac{1}{\mathsf{\tau }}$$

(18)

$${\mathsf{\pi }}_3=\frac{\mathsf{\alpha }\mathrm{h}\mathsf{\tau }}{\mathrm{R}\mathrm{k}}$$

(19)

Since the unknown τ appears in both the monomial π₁ and π₃, and it is interesting that it is only found in one of them, a new monomial can be obtained without this unknown by means of simple mathematical operations between both monomials.

$${\mathsf{\pi }}_4=\frac{\mathrm{h}\mathrm{R}}{\mathrm{k}}$$

(20)

If symmetry had not been assumed, the π₄ monomial would have been ${\mathsf{\pi }}_4=\frac{\mathrm{h}\mathrm{D}}{\mathrm{k}}$, where D is the cylinder diameter, which is the well-known Nusselt number, Nu, which relates the heat transfer coefficient, h, and the material thermal conductivity, k [39]. This same monomial could have been deduced by nondimensionalization from the equivalence between Equations (6) and (8) on the boundary [26].

Finally, the solution for the unknown τ, that is, the time in which the system reaches a steady state, can be deduced by applying the π-theorem, π₁ = Ψ (π_2, π₄).

$$\mathsf{\tau }=\frac{{\mathrm{R}}^2}{\mathsf{\alpha }}\mathsf{\Psi }\left(\frac{{\mathrm{R}}^2}{{\mathrm{H}}^2},\frac{\mathrm{h}\mathrm{R}}{\mathrm{k}}\right)$$

(21)

where Ψ is an unknown function. Analysing Equation (21), and as in the case study, Section 6, the geometric variables, r₁, r₂ and H, will be the same for all materials. The speed at which the temperature increases for each of the vessels will depend on the thermal diffusivity, α, as expected, and on the relationship between the heat transfer coefficient, h, and the material thermal conductivity, k, that is, the Nusselt number. Finally, as previously indicated, the properties of the materials in this study depend on temperature, so this dependence should have been applied to the previous procedure. However, the relationships obtained between the properties of the materials would have been similar, being qualified by their dependence on temperature [21].

5. Material Properties and Model Validation

5.1. Material Properties Depending on Temperature

As mentioned above, we are going to focus, as an example, on the study of the heat transfer that takes place in a water bottle, considering that the value of the density, specific heat, and thermal conductivity of the materials depend on the temperature [10,40,41,42,43]. The materials recently used for water bottles are aluminium (Al319 in this case), polyethylene terephthalate (PET) and polypropylene (PP).

The following equations have been used to determine the density, thermal conductivity, and specific heat of the aluminium 319 (Al319) that will be the material of which our container is made:

$${\mathsf{\rho }}_{\mathrm{A}\mathrm{l}}\left(\mathrm{T}\right)=2668.4118-0.3111\mathrm{T} (\mathrm{kg}/{\mathrm{m}}^3)$$

(22)

$${\mathrm{k}}_{\mathrm{A}\mathrm{l}}\left(\mathrm{T}\right)=76.64+0.2633\mathrm{T}-2·{10}^{-4}{\mathrm{T}}^2 (\mathrm{W}/\mathrm{m}\mathrm{K})$$

(23)

$${\mathrm{c}}_{\mathrm{e}\mathrm{A}\mathrm{l}}\left(\mathrm{T}\right)=747.3+0.2\mathrm{T}+5·{10}^{-4}{\mathrm{T}}^2 (\mathrm{J}/\mathrm{k}\mathrm{g}\mathrm{K})$$

(24)

The properties are represented graphically as follows, Figure 4:

The inside of the bottle will contain water, so it is also necessary to know the above-mentioned properties of water. The following equations have been used to determine the density, thermal conductivity and specific heat:

$${\mathsf{\rho }}_{\mathrm{W}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}\left(\mathrm{T}\right)=-0.0046{\mathrm{T}}^2+2.5216\mathrm{T}+656.4 (\mathrm{kg}/m^3)$$

(25)

$${\mathrm{k}}_{\mathrm{W}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}\left(\mathrm{T}\right)=-1·{10}^{-5}{\mathrm{T}}^2+0.009\mathrm{T}-0.9864 (\mathrm{W}/\mathrm{m}\mathrm{K})$$

(26)

$${\mathrm{c}}_{\mathrm{W}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}\left(\mathrm{T}\right)=-2·{10}^{-4}{\mathrm{T}}^3+0.2026{\mathrm{T}}^2-69.268\mathrm{T}+12010 (\mathrm{J}/\mathrm{k}\mathrm{g}\mathrm{K})$$

(27)

Properties are represented graphically as follows, Figure 5:

To determine the density and specific heat of PET which is the traditional material used in water bottles, the following equations have been used.

$${\mathsf{\rho }}_{\mathrm{P}\mathrm{E}\mathrm{T}}\left(\mathrm{T}\right)=-0.6022\mathrm{T}+1038.2 (kg/m^3)$$

(28)

$${\mathrm{c}}_{\mathrm{e}\mathrm{P}\mathrm{E}\mathrm{T}}\left(\mathrm{T}\right)=0.011{\mathrm{T}}^2-2.8893\mathrm{T}+1045.5 (J/kgK)$$

(29)

The thermal conductivity of PET can be considered constant since, in the temperature ranges in which we are going to work, there are no significant differences in its value. The value of the thermal conductivity will be 0.2976 (W/mK).

Properties are represented graphically as follows, Figure 6:

To determine the density, specific heat and thermal conductivity of PP, the following equations have been used:

$${\mathsf{\rho }}_{\mathrm{P}\mathrm{P}}\left(\mathrm{T}\right)=-0.00305{\mathrm{T}}^2+1.6463\mathrm{T}+625.87 (\mathrm{kg}/{\mathrm{m}}^3)$$

(30)

$${\mathrm{k}}_{\mathrm{P}\mathrm{P}}\left(\mathrm{T}\right)=-\mathrm{0,0016}\mathrm{T}+0.6872 (\mathrm{W}/\mathrm{m}\mathrm{K})$$

(31)

$${\mathrm{c}}_{\mathrm{e}\mathrm{P}\mathrm{P}}\left(\mathrm{T}\right)=0.1417{\mathrm{T}}^2-72.746\mathrm{T}+\mathrm{11,219} (\mathrm{J}/\mathrm{k}\mathrm{g}\mathrm{K})$$

(32)

Properties are represented graphically as follows, Figure 7:

As shown in the figures above, in many cases, the properties show wide variations in their value with temperature.

5.2. Model Validation

To validate both the mathematical model and the network model, a laboratory test was carried out with three cylindrical bottles of similar dimensions made with the materials indicated in the previous subsection,

Table 1. Bottle dimensions.

Material	Al319	PET	PP
Length (cm)	20.5	22.0	23.0
Radio (cm)	3.50	2.75	3.50
Thickness (mm)	3.05	1.80	2.20

. It should be noted that the aluminium bottle was covered with a thin layer of insulating material, which in the network model has been implemented as a resistance. The experiment was carried out in a room that allowed the total absence of light without the influence of any type of radiation and with a constant room temperature during the entire test. Thus, the heat transfer was carried out mainly by conduction and convection. Therefore, the control of the input variables facilitates the validation of the proposed model. To measure the water and environmental temperatures, a thermometer with an error of ±0.1 °C was used, and the measurements were always taken in the central part of each of the cylindrical bottles. In order to compare these results with the proposed model, the properties that influence the problem (density, thermal conductivity, and specific heat) were introduced depending on the temperature at which the material was found. This relationship with temperature was obtained from the references mentioned above.

The experimental procedure consists of introducing water at a temperature higher than 50 °C in each of the bottles and exposing them to cooling conditions, mainly by convection, and without heating due to solar radiation. Temperature measurements were made at the midpoint of the vertical axis of each of the bottles and measured every 15 min until reaching a total time of one hour,

Table 2. Temperature measurements.

	Temperature (°C)
Time (Minutes)	Al319	PET	PP
0	54.1	66.7	55.7
15	53.6	58.7	50.7
30	53.2	52.5	46.7
45	52.5	48.0	43.3
60	52.2	43.6	40.3

. The ambient temperature at the beginning of the tests was 18.7 °C.

As stated above, the experimental tests were simulated by introducing the properties described above for each of the materials. Figure 8, Figure 9 and Figure 10 show the comparison between the results obtained in the simulation and those obtained experimentally. As can be seen, the trend obtained validates the proposed mathematical model, especially in the case of aluminium. In the cases of plastic materials, PP and PET, the deviation is greater due to the heterogeneity in the properties of these materials, as indicated in the previous section. However, in all cases, the evolution trend shown in the simulation is the same as that shown in the experimental test.

6. Results, Case Studies, and Conclusions

Once the mathematical model has been validated, we want to study the behaviour of these materials, used in commonly used water bottles, when exposed to summer temperatures, which in some areas of the Mediterranean can range between 30 and temperatures close to 50 °C in the shade. The foundation of this study is based on the migration of substances from plastic materials to water. Recent studies have determined that at a temperature of 40 °C, spring water in high-density polyethylene (HDPE) bottles contains 180 ng/L of nonylphenol, known for its endocrine-disrupting potential [28]. Other PET studies have found concentrations from 19 to 78 ng/L of this substance in mineral water from PET bottles [44]. Other investigations carried out on plastic water bottles have found migrations of compounds such as dimethyl phthalate (DMP), DEHP, and OP, among others, for PP and PET materials [27].

In this way, in order to compare the behaviour of the materials, it is assumed that the bottles will have the same dimensions, 25 cm in height, 3 cm in radius and 1 mm in thickness. In all cases, the initial temperature of both the water and the material, Al319, PP and PET, will be 20 °C, and they will be exposed to different ambient temperatures,

Table 3. Temperature results of the simulations (°C).

	Temperature at Reference Point (°C)
Room Temperature (°C)	Al319	PET	PP
30	23.72	27.77	27.58
40	27.43	35.59	35.14
50	31.13	43.48	42.69
60	34.72	51.43	50.22

. The heating of the material occurs mainly by convection since we have assumed that they are isolated from solar radiation. Thus, Figure 11, Figure 12 and Figure 13 show the temperature distribution for one hour for each of the materials and ambient temperatures. The black line on the right side of the image represents the boundary between the water and the bottle material. This point will be the reference point for comparisons between the materials studied, Table 3. In addition, the results for an ambient temperature of 60 °C have been added. As can be seen, the material that increases its temperature the least is Al319, presenting temperature differences of up to 12 °C with plastic materials for an ambient temperature of 50 °C, possibly because of the thin layer of insulation on the bottle. Plastic materials, PET and PP, present similar temperature values, being slightly lower in the case of PP. In both cases, we find temperatures close to ambient. It should be noted that the results obtained in this test depend on the materials used, and other aluminium alloys or variations in the composition of plastic materials could significantly change the result obtained. It is important to highlight that the behavior of the bottles validates the results obtained through the nondimensionalization technique, Equation (21). Finally, in all the distribution figures, the radial distribution of the problem in temperatures is clearly appreciated.

In conclusion, in this paper, a mathematical model capable of simulating cylindrical containers containing fluids whose properties are dependent on temperature has been fully developed. The simulation has been carried out using the Network Simulation Method, where the construction of the network model has been presented in detail. In addition, this has been validated with experimental tests carried out on three cylindrical containers made of aluminium, PET, and PP.

Moreover, the properties of materials whose value depends on temperature have been presented, highlighting the difficulty, in some cases, of obtaining experimental values as a function of temperature in the bibliography since it is common to take a constant value of these properties [10]. However, nowadays, more and more works are showing the dependence of these properties on temperature, making it necessary to develop mathematical models that include this dependence, thus improving the reliability of the simulations [10,11,12,13,14,15,16].

On the other hand, the nondimensionalization technique has also been applied, which allows us to obtain the weight of each of the variables on the problem, determining that the cooling or heating speed depends on the thermal diffusivity, as expected, and of the relationship between heat transfer coefficient and the material thermal conductivity, that is, the Nusselt number.

Finally, by way of illustration, a case has been studied where aluminium, PET, and PP mineral water bottles are exposed to summer environmental conditions, where temperatures close to 50 °C are reached. It should be noted that the best behaviour is presented by aluminium with a temperature difference of up to 12 °C with plastic materials, possibly because of the thin layer of insulation on the bottle. As a conclusion of the case study, the temperatures reached by the water in the summer season could imply migrations of different substances from plastic materials, PET and PP.