An Assessment of the Dynamic Behaviors of a Counterflow Heat Exchanger Under Variable Thermophysical Properties of Thermal Fluids

Heat exchangers (HEs) find applications across various sectors and systems, including food industries, refinery processes, renewable energy systems, waste heat recovery, automotive engines, and medical applications. The size and configuration of HEs vary based on specific requirements and intended applications. Various heat exchanger configurations, including crossflow, shell-and-tube, and plate heat exchangers, facilitate the transfer of thermal energy between different media. Through them, plate heat exchangers are extensively used in many industrial areas due to their positive contribution to processes, the economy, and the environment [1, 2]. Heat exchangers involve three heat transfer mechanisms: conduction, convection, and radiation. In general, conduction and convection heat transfer mechanisms generally dominate the heat transfer processes due to the assumptions of adiabatic processes and well insulation from the environment [3]. Researchers work to improve heat transfer by altering the plate geometry in an effort to make them more efficient and compact [4]. Therefore, before constructing heat exchangers, the interaction of fluids and the wall separating the fluids should be studied [5]. Especially in transient operating conditions, the change in the temperature distribution inside the heat exchanger and determining the time taken for the transition to a stable state are significant for the healthy operation of systems [6, 7]. For this reason, intensive research continues using different methods for different types of heat exchangers to determine the behavior of heat exchangers under dynamic conditions. For this purpose, various algorithms and models were developed using different software and simulation tools such as C#, MATLAB/SIMULINK, ANSYS-FLUENT, COMSOL, etc.

Numerical calculation methods provide advantages such as rapid and easy applicability, flexibility, high accuracy, and low hardware and time requirements. Owing to these benefits, MATLAB/SIMULINK provides a graphical programming environment for dynamic and complex models and simulations and is employed in many engineering and scientific studies [8]. For instance, Bobič et al. [9] proposed a sophisticated SIMULINK model specialized for corrugated counterflow plate heat exchangers to observe the temperature changes concerning time under transient conditions. They performed experimental studies to validate the model. They followed a good agreement between practical and theoretical studies. Da Silva and Fernandes [10] conducted a study including thermodynamic modeling of a PV/T solar system parametric simulation using MATLAB/SIMULINK. They developed the thermodynamic model of the PV/T system using a differential heat transfer equation for the solar collector and storage tank separately. Results showed that the total energy efficiency of the PV/T was found to be 24%, 15% thermal, and 9% electricity. Ilis et al. [11] studied the construction and working parameters of an adsorption chiller using MATLAB/SIMULINK simulation. The flow rates of the flows in an adsorber and chilled water, the heat transfer surface area of the evaporator, and the material of the chiller were considered input parameters that affect the heat transfer. They revealed that the most promising results were observed by aluminum material if the water flow rates supported the performance. Bologa et al. [12] performed a parametric study using MATLAB/SIMULINK to investigate the heat transfer and storage performance of the Intermediate Heat Transport and Storage System. They considered tank geometry and temperature of molten salt as input data that affects temperature decrease in the system tanks due to the service time. As a result, they revealed that increasing tank height with constant molten salt mass decreases the duration required for the molten salt to reach critical working temperature.

Thermal systems experiencing heat exchange with one or two fluids undergoing phase change are an example of dynamic heat exchangers. A wide range of thermal systems use phase change materials (PCM) as sources or sinks of heat, like Paraffin wax and Molten salt. Al-Kayiem and Alhamdo [13] performed experimental investigations and developed in-house code to simulate Paraffin wax melting and solidification as TES material in a packed encapsulated heat exchanger. The developed code is a transient, one-dimensional numerical model. Deng et al. [14] proposed a SIMULINK model for a plate heat exchanger consisting of a phase change material (PCM) interface based on the concept of printed surface heat exchangers and tested the dynamic responses. They revealed that the usage of PCM prevented efficiency decrease under temperature disturbances. Trafczynski et al. [15] proposed a SIMULINK model for PID-controlled shell and tube heat exchangers (STHE) to discover the effects of fouling on dynamic behaviors. The heat exchanger was divided into sub-cells in the study, and a mathematical model was created using the lumped model approach. Results showed that according to the thermal resistance of fouling, tuning parameters have to change.

It has been realized that in previous studies, SIMULINK models are intensively employed to investigate the complex behaviors of various heat exchangers under instant temperature disturbances. However, the studies were conducted with the constant thermophysical properties assumptions for fluid flows and heat exchanger materials. Many fluid flows possess changes in the thermophysical properties, in particular, due to temperature change, which is a point of research gap.

In this study, thermal interactions in transient operating conditions of fluids and the wall at crossflow heat exchangers (CFHEs) operated with ethylene glycol and engine oil are determined by the MATLAB/SIMULINK model considering the variable thermophysical properties to observe the effects of the temperature difference and the errors due to the constant temperature assumption with disturbances in the inlet temperatures. Differential equations were written to develop the thermodynamic model of the counterflow plate heat exchanger that was divided into n cells imaginarily. Then, equations defining heat and mass transfer were considered for two-dimensional heat transfer: between hot flow, cold flow, and heat exchangers' walls, regarding the variation of thermophysical properties of flows and heat exchanger materials by temperature. Hence, the differentiation of temperature distributions of cells in heat exchangers was instantly observed under transient operating conditions. Finally, fluid and wall temperatures are presented in graphs, and how instantaneous changes affect the heat exchanger's temperature distribution.

In this study, the two-dimensional heat transfer phenomenon of a CFHE made of AISI 304L stainless steel operated by engine oil on the hot flow side and ethylene glycol as cooling liquid on the cold flow side is considered. Firstly, a thermodynamic model was prepared for the CFHE, and the MATLAB Simulink tool was used to solve the differential heat transfer equations in the thermodynamic model. Thermophysical properties of AISI 304L material were received from Kim [16]. Engine oil and cooling liquid ethylene glycol thermophysical properties were acquired from Incropera et al. [17].

First, the heat exchanger was divided into n cells and examined as inlet cells, middle cells, and outlet cells. Thus, the differential equations defining the heat and mass transfer were created separately for the first cell, which is the input cell, the cells between the first cell and the nth cell, and finally, the n'th cell, considering the finite volume method (FHV) described in Aragón and Duarte [18] and Salazar-Herran et al. [19]. Thus, the heat transfer model in the current study consists of five environments, namely the lower wall, cold fluid, middle wall, hot fluid, and upper wall. In addition, separate differential equations are created for three cells: the first cell is any i'th cell between the first and the n'th cells. Thus, for the CFHE, a total of fifteen equations for n cells were solved with the help of MATLAB SIMULINK. The general view of the heat exchanger subject to this study, the flows' directions and locations, and the order of cell numbers are demonstrated in Figure 1.

In the thermodynamic model, potential and kinetic energy changes, the pressure drop between cells, heat loss to the environment, radiation heat transfer, and phase changes in fluids are ignored. In addition, temperatures are evenly distributed in each cell, and fluids fill the cells. Thermophysical properties such as heat conductivity, specific heat, dynamic viscosity, density, and Prandtl number are simultaneously calculated for differing temperatures across the heat exchanger using the thermophysical properties versus temperature variation supplied by Aragón and Duarte [18] inserted into the model.

In addition, heat conductivity, specific heat, and density values are calculated using the correlations derived as a function of temperature as given in Eqs. (1)-(3) for considered heat exchanger material AISI 304L stainless steel [16]. Hence, temperature change for all thermophysical properties is considered in this simulation model for fluids and walls. In Figure 1, L, w, and h indicate width, the height of the channels, respectively, and t is the thickness of the walls.

Differential equations consisted of heat and mass transfer under transient conditions for the lower wall, cold flow, inner wall, hot flow, and upper wall at the first cell, given in Eqs. (4)-(8) from the lower wall to the upper wall.

For i=1,

Differential equations for heat and mass transfer under transient conditions for lower wall, cold flow, inner wall, hot flow, and upper wall at i'th cell are given in Eqs. (9)-(13) from lower wall to upper wall, respectively.

For i>1 and i

Differential equations consisted of heat and mass transfer under transient conditions for lower walls, cold flow, inner wall, hot flow, and upper low at n'th cell in Eqs. (14)-(18) from the lower wall to the upper wall.

For i=n,

where, m (kg), ṁ (kg/s), c (kJ/kg‧K), h (W/m²‧K), k (W/m‧K), A (m²), T (℃), L (m), and n indicate mass, mass flow rate, specific heat, heat convection coefficient, heat conduction coefficient, heat transfer surface area, temperature, length of the heat exchanger and the number of divisions of heat exchangers respectively. The subscripts w, l, c, css, n, and t indicate the wall, lower, cold, cross-section, inner, and top, respectively. The thermal interactions between the 1st and nth cells are shown in Figure 2, and according to Eqs. (19)-(23).

In the second step, the heat transfer coefficients of hot flow, h_h, and cold flow, h_c, used in the differential equations were calculated according to the parallel flow principles on the flat counterflow as given in Eqs. (24)-(27), taking into account the change of the thermophysical properties of the fluids with temperature [17]. Finally, Reynolds' numbers for hot flow and cold flow were calculated by Eq. (24).

Average Nusselt numbers over the entire counterflow for laminar and turbulent flow were used in this study. Under laminar flow conditions (Re_L <5×10⁵), the Nusselt number was calculated using Eq. (25) [20].

Under turbulent flow conditions (5×105≤Re_L≤107 and 0.6≤Pr ≤60) Nusselt number was calculated using Eq. (26) [20].

Heat convection coefficients for both fluids are calculated using Eq. (27) [20].

Geometric dimensions in the current model for the CFHE, such as wall thickness (t), counterflow length (L), channel height (h) and width (w), and thermophysical properties, such as heat conduction coefficient (k), dynamic viscosity (μ), specific heat (Cp) mass flowrates of the flows (ṁ) of the flows are determined as input parameter for the Simulink model. Temperatures of the hot flow, Th, cold flow, Tc, and walls, Tw, are defined as output parameters.

The operating principles of the model in Simulink are summarized in Figure 3. The initial temperatures of the fluids and the heat exchanger material, the mass flow rates of the fluids, and the geometric information of the channels through which the fluids pass were determined as inlet conditions. Reynolds and Nusselt's numbers are calculated with the help of MATLAB by determining the thermophysical properties of the fluids and the heat exchanger material according to the inlet conditions information. The heat transfer coefficients on both sides are calculated using Reynolds and Nusselt's numbers. The results obtained are transferred to the SIMULINK block diagram given in Figure 4, and the cell outlet temperatures of the fluids are calculated. These results establish the entry conditions for the next cell. The thermophysical properties of the fluids are recalculated according to the new entry conditions and the calculations are repeated until the last cell.

Figure 4 shows the block diagram of the Simulink model in which the cell outlet temperatures are calculated. Inlet conditions of the lower inner and upper walls were selected at 10℃, and the flow temperatures were selected at 20℃ and 80℃ for hot and cold flows. 1/s stands for the derivative of the dt/dt of the initial conditions. The heat transfer coefficients calculated in MATLAB are substituted in the differential equations created for the cell model, and they are solved with the fourth-order Runge Kutta integration method, which is an excellent differential equation solver, as given in Eqs. (28)-(32) [21]. Finally, the cell exit temperatures are calculated.

The effective slope used is a weighted average of the slopes at the four points $\left(x_n, y_n\right)$, $\left(x_n+\frac{h}{2}, y_n+\frac{1}{2} k_1\right)$, $\left(x_n+\frac{h}{2}, y_n+\frac{1}{2} k_2\right)$, and $\left(x_{n+1}+\frac{h}{2}, y_n+k_3\right)$ in the x-y plane, an average because the sum of coefficients 1/6, 2/6, 2/6, 1/6 that multiply the k's is equal to 1.

The following equation was used to calculate the difference between the thermophysical properties that change depending on the temperature with the assumption of constant temperature [22].

Time-dependent MATLAB-SIMULINK simulations were carried out for the CFHE by using constant and temperature-varying thermophysical properties. In addition, a short-term disturbance was applied for variable thermophysical properties, and the responses of the CFHE were observed.

Figure 5 shows the temperature distribution of the CFHE versus cell numbers. According to Figure 5, the ethylene glycol entered the heat exchanger from 1st cell at 10℃, leaving the heat exchanger from the 1000'th cell in different temperatures, such as 64℃ and 72℃ for constant and variable properties, respectively. Similarly, the engine oil enters the heat exchanger from the 1000th cell at 80℃, leaving the heat exchanger from the first cell at different temperatures, such as 22℃ and 37℃, for constant and variable properties. Temperature differences in engine oil for variable and constant properties are observed at 43℃ and 58℃, respectively. Thus, the simulation using fixed properties gave an erroneous result of 15℃ at the hot flow exit temperature. Temperature differences in ethylene glycol for variable and constant properties are observed at 72℃ and 64℃, respectively. Thus, the simulation using fixed properties gave an erroneous result of 8℃ at the hot flow exit temperature.

Effectiveness values for variable properties and constant properties are calculated at 0.614 and 0.828, respectively. Although the temperature difference is relatively small, the ethylene glycol temperature change across the heat exchanger is very different. This difference is due to the high variation in thermophysical properties of ethylene glycol, even at minor temperature differences. For example, while the specific heat of water is 4.198kJ/kg‧K at 280 K, there is only a maximum 0.38% change at 370 K with 4.214kJ/kg‧K. However, the specific heat of ethylene glycol and engine oil varies by 17.43% and 20.74%, respectively, in the same temperature range. This change is clearly seen in the temperature distribution that the flows show across the heat exchanger. Accordingly, it indicates that constant properties should not be assumed in the flows of fluids whose thermophysical properties vary significantly according to temperature, such as ethylene glycol and engine oil.

Figure 6 demonstrates the temperature distribution of the heat exchanger's upper, inner, and lower walls. The upper and lower wall temperature distributions versus cell numbers are observed to be identical with the hot flow and cold flow temperatures at the end of 20 s simulation time. For instance, inner wall temperatures for variable and constant properties are observed at 18℃ and 15℃ for the 1st cell and 75℃ and 70℃ for the 1000th cell.

Differences between the variable properties and constant properties assumptions versus temperatures for the ethylene glycol and engine oil are shown in Figure 7. Errors in temperatures of ethylene glycol are increasing sharply from zero to 25% up to the 180th cell, and it shows a slight decrease to 22% up to the 400th cell and repeating the increase to 25% up to the 400th cell. After, it offers a rapid decline to 12% at the 1000th cell. Inner wall temperature error follows a relative path to the cold flow errors, and it finished with 6% in the 1000th cell.

However, temperature errors in hot flow show continuously increasing values from the 1000th cell to the 1st cell. Next, errors show a rapid increase from zero to 24% up to the 600th cell. Then, the growth of the error becomes slower and reaches 27.5% up to the 350th cell. Finally, however, it shows a sharp increase from 27.5% to 42.5% up to 1st cell.

The error rates vary between 20-30% for hot flow inner wall and cold flow, while cell numbers are between 200-750. Finally, the constant properties assumption caused a 42.5% error for hot flow outlet temperature and 12% for cold flow outlet temperature. Maximum errors are observed for hot flow at the outlet (1st cell) with 42.5%, cold flow at the 600th cell with 25%, and the inner wall at the 75th cell with 25%.

Figure 8 shows the application duration, time, and level of temperature disturbance for hot flow. For example, temperature disturbance is applied for the inlet temperature at the 4th second for 2 seconds with a 60℃ increase, and after the 6th second inlet temperature returned to 80℃.

Response of the heat exchanger is observed with variable properties in Figure 9. The temperature distribution is generated depending on the cell number and time domain up to 1000th cells and 20 s, respectively. According to the 3-D surface response graphs in Figure 9 (a1), the temperature distribution for cold flow increases by increasing duration and cell number in the first five seconds and 200 cells rapidly. However, it is seen that the rate of growth in cold flow temperature decreases in the following seconds. When the cold flow meets the hot flow, it is seen in Figure 9 (a1) that the sudden temperature increase is seen in the 2nd second. Therefore, in the next seconds, the temperature increase decreases as the fluid temperatures approach each other. Figure 9 (a2) shows the reaction of the cold flow when a sudden temperature change of 60℃ is applied to the hot flow. According to this, it is seen that the temperature in the 1000th cell, where the hot flow enters, increases rapidly from the 4th second to 82℃ after 2 seconds. Then the temperature decreases slowly to 74℃ at the end of the 20th second. According to Figure 9 (b1), the hot flow inlet temperature drops suddenly from 80℃ to 62℃ in the first 3 seconds. This decrease is the high heat flux due to the wall temperature at 25℃ in initial conditions. In the following seconds, the temperature of the hot flow increases up to 78℃ with the increasing wall temperature, as seen in Figure 9 (c2), and decreases to 37℃ in the 1st cell at the end of the 20th second by transferring heat to the cold flow.

According to Figure 9 (b2), the sudden temperature increase in the 4th second caused a rise in the first cell from 78℃ to 128℃ as a result of the heat flux to the wall and cold flow. At the end of the 5th second, the temperature increased and reached 130℃. Then, the hot flow temperature dropped to 42℃ in the 1st cell after 20℃. In Figure 9 (c1), the interior wall temperature shows a similar variation characteristic to the cold flow temperature. Accordingly, the internal wall temperature increased from 25℃ to 38℃ at the end of the first second.

When examining existing literature, it has been observed that this method is commonly employed in heat exchange processes involving fluids with varying thermophysical properties [22]. It is particularly relevant in scenarios with significant temperature differences and in heat exchangers handling cryogenic fluids [23, 24]. Additionally, accurate calculations for low-temperature heat recovery systems, such as CO2 cycles, must account for these variable features. The study highlights that even in low-temperature ranges, fluid thermophysical properties can significantly impact heat transfer calculations [25].

In this study, a dynamic simulation for counterflow heat exchangers (CFHEs) is performed using a MATLAB/Simulink model. The temperature distribution of hot flow, cold flow, bottom, inner, and top wall in counterflow is simulated under varying conditions due to the distortions in inlet temperatures. It was revealed that the assumption of constant thermophysical properties in heat exchanger analyses is not suitable for some fluids such as ethylene glycol and motor oil. As a result of time dependent simulations made according to constant and variable thermophysical properties, the error ratio was observed that it was approximately 43% in hot engine oil and approximately 23% in ethylene glycol. In addition, the temperature distributions of heat exchangers with complex geometries were obtained according to variable inlet conditions and time-dependent at constant inlet conditions.