*3.1. General Description*

In Section 1, it was mentioned that the tested PCM is dedicated to working in a special storage heat exchanger/tank [38,39]. It was decided that the PCM, due to the design of the heat exchanger, should be studied in a cylindrical system for the conditions most similar to those in which it will be operated later. For this purpose, a wind tunnel, as presented in Figure 4, was used that met the operational and functional parameters. The obtained range of the Reynolds numbers (15,260–52,767) is fully turbulent and the heat transfer conditions (temperature and the Nusselt numbers) are similar to the considered storage tank.

**Figure 4.** General view of the wind tunnel with installed PCM module.

The experimental set-up consists of four basic sections, operated in an open circuit suction mode. The first section is related to the air inlet. It is composed of a lamellar exchanger with a cross-section of 0.75 m × 0.74 mm in which it was possible to control the temperature in the range of −15 ◦C to 35 ◦C while using a chiller with a maximum cooling capacity of 4 kW. Subsequently, a stabilisation flow section was installed behind the lamellar exchanger. The stabilisation section is made of three net layers with a mesh size of approximately 1.5 mm × 1.5 mm. Such a system ensured a stable flow, without recirculation zones while also eliminating potential turbulence behind the lamellar exchanger. The next section is the so-called the Witoszy ´nski nozzle, whose task is to shape a flat velocity profile in the measuring section. This special velocity profile allowed for homogeneous thermal-flow

parameters along the entire length of the tested module to be obtained. The measuring section had a square cross section with dimensions of 0.25 m × 0.25 m and a length of 0.52 m. In order to eliminate the impact of the installed fan, between the measuring section and the fan a 0.9 m long channel with an identical cross-section as the measuring section was installed. The purpose of this channel was to stabilise the fan working conditions and eliminate any possible flow disorders that are caused by its operation. The fan that is installed in the experimental set-up provides the opportunity to achieve three mean values of velocity in the measuring section equal to 0.92, 2.27, and 3.18 m/s, respectively, which correspond to the Reynolds numbers 15,260, 37,688, and 52,767, based on the channel height of the test section. In order to achieve different mean value of air velocities, the fan power was regulated by a capacitor system. Generally, the maximum volume air flow rate was 1000 m3/h, which was generated for the maximum electrical power of the fan engine, equal to 150 W.

The parameters of inlet air, such as ambient temperature, pressure, and humidity, were measured just before the lamellar exchanger by a high-precision digital temperature, humidity, and airflow meter Testo 480 with Testo Robust Humidity Probe.

## *3.2. Tested Module*

The PCM was experimentally tested in a special cylindrical module made of copper, as presented in Figure 5. This module was characterised by a fixed height and a wall thickness of 250 mm and 1 mm, respectively. Three different outer diameters of 15, 22, and 28 mm were studied. At the ends of the module, two pipe caps that were made of the same material with a height of 15 mm were mounted. Additionally, three T-type thermocouples have been installed in situ in the axis of the module, at half height, and one-third above and below, as is shown in Figure 5. The module was then placed and tested in the measuring section, as depicted in Figures 4 and 6.

**Figure 5.** Tested module of a PCM.

The PCM mass contribution in the module varied from 25.23% to 37.82%, (see Table 2). The theoretical heat transfer from the air to the modules were estimated and presented in Table 2, according to the producer information, as described in Table 1 and assumptions concerning paraffin and copper properties used during the numerical calculations, as summarised in Table 4. For paraffin, it varies from 2.87 kJ to 10.19 kJ for sensible heat and from 2.46 kJ to 8.74 kJ for latent heat. For copper the heat transfer varies from 1.62 kJ to 3.71 kJ. The heat storage in copper is 23.31% for *d* = 15 mm, 20.78% for *d* = 22 mm and 16.37% for *d* = 28 mm of the total system, respectively. Those estimations show that PCM phase change plays an important role during the heating of the module that is filled with the PCM.


#### *3.3. Preparation of the Module to the Test*

The module has been specially prepared in order to test the thermal behaviour of the module filled with the PCM for various thermal-flow conditions. First, before filling with the PCM, the module was degreased with propyl alcohol and thoroughly dried. Next, three T-type thermocouples with an accuracy of ±0.5 ◦C were placed in the prepared arrangemen<sup>t</sup> in the configuration that is shown in Figure 5. Subsequently, the module was filled with the phase change material. Before testing, the module was placed in a freezer for about 4 h to cool to a temperature of about −18 ◦C. The freezing process was completed when the material obtained the same temperature in the entire volume, which was monitored by thermocouples that were placed inside the module. The module prepared in this way was then mounted in the wind tunnel measuring section. Figure 6 presents the measuring section with the installed PCM module.

**Figure 6.** PCM module placed in the wind tunnel.

#### *3.4. Measurement and Control Systems*

The experimental set-up was equipped with two types of measurement system: one was based on Testo devices and the second one on the National Instruments module 4-Slot USB CompactDAQ Chassis (cDAQ 9174) working under the LabView environment. The high-precision digital temperature, humidity, and airflow meter Testo 480 had a buildin absolute pressure meter with an accuracy of ±3 hPa and allowed for up to three probes to be connected. The first Testo device was a Robust Humidity Probe, which was used to measure humidity and temperature of the suction air with an accuracy of temperature ±0.5 ◦C (from −20 ◦C to 0 ◦C) and ±0.4 ◦C (from 0.1 ◦C to 50 ◦C), relative humidity: ±2% RH (from 2.1 to 98% RH). The second installed device was a Thermal Flow Velocity Probe, Testo 480, which was used to measure the velocity profile. The device is characterised by the accuracy of measured velocity profile amounting to ±0.03 m/s + 5% of the measured velocity.

A special thermocouples net system was installed in order to measure inlet and outlet air temperatures at the measured section. At the inlet to the section, eight T-type thermocouples were installed and used to determine the average inlet temperature. At the outlet, one T-type thermocouples type of temperature sensor was mounted. The accuracy of installed T-type thermocouples was ±0.5 ◦C.

The combination of the National Instruments equipment and LabView software allowed for the measurement of temperature with a high resolution and frequency. During the experiments, it was sufficient that the measurements were carried out with a frequency of 0.2 Hz (every 5 s).

#### *3.5. Velocity Profiles at the Inlet and Outlet of the Measurement Section*

A special construction of the inlet section with the installed Witoszy ´nski nozzle contributes to the formation of quasi-uniform thermal and velocity profiles in the measurement section, as mentioned in Section 3.

Before the measurement campaign, the velocity profiles at the inlet and outlet were determined. The measurements of the velocity profile were carried out at a distance of 20 mm from the inlet and outlet edges of the measuring section. The velocity profiles were determined at measuring points 33, 73, and 105 mm from the channel's symmetry axis.

In total, seven holes have been made in the top of the measuring section cover for measurement purposes. The thermal flow velocity probe that was connected to the Testo 480 module has been inserted into these holes, mounted on a special measuring instrument. A measuring device, similar to a caliper allowed for precise measurement of the velocity profile along with the channel height with an accuracy of ±0.1 mm. In total, one profile was determined while using 19 measurements, concentrated in the central part of the channel. The results of the obtained profiles at the inlet and outlet of the measuring section for different powers of the fan are presented in Figure 7, where a larger number represents higher air velocity.

The measurements that are presented in Figure 7 showed that, for all the tested ranges of fan power, the velocity profiles are characterised by a flat profile. The obtained flat velocity profile is a consequence of the use of a specially shaped inlet contraction section designed according to the Witoszynski curve equation [41]. The analysis points out that the differences between the average values determined do not exceed 7.4% of the local value of velocity. Table 3 summarises the test results. In addition, for further analysis, each analysed value of average speed: 0.0, 0.92, 2.18, and 3.18 m/s was defined as Case 0, Case I, Case II, and Case III, respectively.

**Table 3.** The summarised results that were obtained from profile measurements.


**Figure 7.** Velocity profiles of cross section at different power.

#### **4. Numerical Modeling of the Analyzed System with PCM**

*4.1. Process Governing Equations*

The thermal-flow behaviour in the PCM and flowing air can be described by a standard set of equations [42,43] and for the three-dimensional case, comes down to solutions of Partial Differential Equations, (PDEs), which include a mass, a momentum, and an energy equation [44]. In this context, for solving the set of equations, one of the most frequently used numerical methods, called the Finite Volume Method (FVM), [45,46] was applied.

In order to solve transport equations of mass, momentum and energy, the commercial CFD tool-Ansys Fluent has been utilised [47]. To take the phase change of the PCM into account, an enthalpy-porosity technique [48] has been used, which is already implemented in Ansys Fluent software module. Using this technique, the interface is not explicitly tracked, but the liquid fraction for each domain cell is solved in each iteration. The cells that contain a liquid fraction *β* between 0 and 1 create an interface region, called the mushy zone. This region is treated as a porous medium with porosity ranging from 1 for a liquid to 0 for the solid. The porosity is equal to the liquid fraction.

The liquid fraction is calculated knowing the balance of the enthalpy *H*. The enthalpy of a material is computed as a sum of sensible enthalpy *h* and the latent heat content Δ*H*

$$H = h + \Delta H \tag{1}$$

where sensible enthalpy *h* is defined as follows

$$h = h\_{ref} + \int\_{T\_{ref}}^{T} c\_p \mathbf{d}T \tag{2}$$

*href* and *Tref* are reference to enthalpy and temperature, respectively, and *cp* is specific heat at constant pressure. Regarding the liquid fraction, it is defined in the following manner:

$$\beta(T) = \begin{cases} 0 & \text{if } T < T\_{\text{soliding}} \\ \frac{T - T\_{\text{soliding}}}{T\_{\text{liquid}} - T\_{\text{soliding}}} & \text{if } T\_{\text{solulus}} < T < T\_{\text{liquid}} \\ 1 & \text{if } T > T\_{\text{liquid}} \end{cases} \tag{3}$$

Knowing the latent heat L of the material used, the latent heat content in Equation (1) is calculated as

$$
\Delta H = \beta \text{ L} \tag{4}
$$

where Δ*H* varies between 0 for a solid and *L* for a liquid. Based on the liquid fraction distribution, the enthalpy is calculated and then inserted into the energy equation, where the temperature distribution is solved via an iterative manner.

#### *4.2. Geometry, Boundary and Initial Conditions and Numerical Schemes*

The geometry used during numerical calculation is based on the experiments that are described in Section 3. For the purpose of simplifying the calculations, only the measuring section is modelled. Figure 8 presents a scheme of the studied case. During the simulation, to the rectangular duct, made of plexiglass, hot air of a temperature 24 ◦C inflows according to the experimentally obtained uniform velocity changing in a range of 0.92 ÷ 3.18 m/s. The duct is 0.7 m long and the inlet surface is a square, with a side length of 0.25 m. The wall thickness of the duct is 3 mm. In the middle of the duct, the copper cylinder with an outer diameter of 15 ÷ 28 mm and 250 mm long is set. To close the module, pipe caps that were made of copper, 15 mm long and 15 ÷ 28 mm in diameter, were mounted at the ends of the test element. The cylinder and bases are hollowed and the wall thickness is 1 mm. The inside of the cylinder is filled with the tested PCM. The air outflow has an atmospherical pressure of 1 bar.

For the air, plexiglass and copper materials constant thermophysical properties are assumed. For the real PCM, physical properties, such as heat capacity or thermal conductivity, vary with temperature [49]. Thermal conductivity *k* is higher in the solid than in the liquid phase. In the case of the PCM, the density, thermal conductivity, and dynamic viscosity are regarded as temperature-dependent and calculated via Equations (5)–(7), respectively. Table 4 describes the thermophysical properties of all materials that are used during the numerical simulations.

**Figure 8.** Three-dimensional geometrical model used during numerical calculation.

Density of the PCM:

$$\begin{aligned} \text{for } T < T\_{\text{solidus}} &\Rightarrow \rho = 910 \frac{\text{kg}}{\text{m}^3} \\ \text{for } T \ge T\_{\text{solidus}} &\Rightarrow \rho = 790 \frac{\text{kg}}{\text{m}^3} \end{aligned} \tag{5}$$

Thermal conductivity of the PCM:

$$\begin{aligned} for \ T &< T\_{liquidus} \Rightarrow k = 0.5 \,\frac{\text{W}}{\text{m} \,\text{K}}\\ for \ T &\ge T\_{liquidus} \Rightarrow k = 0.4 \,\frac{\text{W}}{\text{m} \,\text{K}} \end{aligned} \tag{6}$$

Dynamic viscosity of the PCM [50]:

$$\eta = 0.001 \exp\left(-4.25 + \frac{1700}{T}\right) \text{Pa s} \tag{7}$$

For all analysed cases, the axis of the cylinder is parallel to the gravity vector. The ratio of the module height to the inside diameter is large and it varies from 8.93 to 16.67, in the PCM volume, natural convection current occurs, which has an important impact on heat transfer and fluid flow behaviour in PCM. To take this phenomenon into account, the Boussinesq approximations that are shown in Equation (8) have been made in the PCM domain using UDF procedure. The following term has been added to the y-momentum equation for cells with liquid fraction higher than 0 within PCM domain:

$$S\_{\mathcal{V}} = -\rho\_{\text{liquidus}} \lg \beta (T - T\_{\text{solidus}}) \tag{8}$$

where: the thermal expansion coefficient of the paraffin equalled *β* = 0.000778 1/K has been included in the numerical calculation procedures [51].

**Table 4.** Thermophysical properties of materials used during calculations.


The geometrical model that is presented in Figure 8 was discretized in space with the use of Ansys Meshing software. A tetragonal mesh was used within the air and PCM space. Solid elements were discretized while using hexahedral elements. Additionally, on the outer and inner surfaces of the cylinder 12 inflation layers have been applied. The thickness of the first element adjacent to the wall was set to 0.1 mm. In order to ensure energy transfer between solid and fluid regions, the applied meshes were non-conformal and were connected with the use of coupled wall boundary condition [47]. Figure 9 presents the employed mesh with details.

**Figure 9.** Mesh with details used during numerical calculations.

Table 5 summarises the complex boundary conditions applied during the simulations.


**Table 5.** Boundary conditions employed during calculations.

Air and PCM were both modelled as incompressible fluids. In order to account for turbulence, in the air region a standard *k* − *ε* turbulence model with enhanced wall treatment was employed. In the PCM region, laminar flow is assumed. The boundary conditions used in the model are depicted in Figure 8. A no-slip boundary condition was utilised on all walls. The initial temperature of the tunnel walls and the air was set to 24 ◦C, while the PCM, cylinder, and lower and upper bases to 2 ◦C.

Calculations have been performed in a transient mode with a time-step of 0.01 s with the first-order implicit treatment. The total time of the simulation was 1440 s. The SIMPLE algorithm with the second-order scheme was used for pressure-velocity coupling. Transport equations of momentum, turbulent kinetic energy, turbulent dissipation, and energy were discretized while using the second order upwind scheme. Gradients were calculated with the use of a least-square cell-based scheme [47]. The solution in each time-step was considered as converged if the residuals were less than <sup>10</sup>−6; however, a maximum of 30 iterations per time-step was done. For energy, liquid fraction, and turbulence equations the under-relaxation of 0.5 was applied. For other equations, the default values were set. The calculations were conducted in a parallel mode on a cluster with the use of 16 Intel Xeon E5-2670 v3 2.3 GHz (Haswell) processors. The time of the calculations was approximately two weeks.
