1. Introduction
PHEs are used in many areas of industry. They are commonly used in heating; refrigeration; district heating; and various industries, such as the food processing, brewing, pharmaceuticals, oil, gas, and chemical industries. These heat exchangers are composed of a number of plates with a thickness of approximately 0.3 (Brazed PHE) or 0.4 mm (Plate-and-Frame PHE) [
1]. Depending on the area of use, heat exchangers differ in the way the plates are connected (brazed, welded, gasketed, and semi-welded) and the material of the plates (stainless steel or more sophisticated alloys, e.g., alloy 254 SMO, alloy C-276, titanium, palladium-stabilised titanium, Incoloy, and Hastelloy) [
1]. Their popularity is due to their small size, but, at the same time, high thermal efficiency. For example, compared with shell-and-tube heat exchangers, the total heat transfer coefficient of the PHE is about twice as high [
2]. The large value of the overall heat transfer coefficient is caused by the embossing on the plates, due to which the fluid flowing between the plates is disturbed, increasing the intensity of heat transfer. The embossing on the plates also increases their stiffness, which allows PHEs to operate at high pressures. The high heat output of PHEs is due not only to the higher values of overall heat transfer coefficients, but also to the fact that the embossing on the plates increases the surface area of heat transfer of these apparatuses.
A significant problem in the operation of PHEs is the efficient assessment of the fouling condition of the apparatus. They should be subjected to regular cleaning; otherwise, fouling will cause a decrease in the heat output of the heat exchanger, as well as an increase in the pressure drop of fluid. Contaminants are deposited in heat exchangers regardless of the fluid physical state [
3,
4,
5]. The literature indicates that in some industries (e.g., the dairy industry [
5,
6,
7,
8]), the occurrence of contaminants in this type of apparatus is a serious problem.
Although fouling in heat exchangers is a problem, as mentioned earlier, in the case of PHEs, due to disturbed fluid flow, there is partial self-cleaning of the apparatuses. This means that PHEs do not foul as quickly as, for example, shell-and-tube heat exchangers. In addition, the deposition of fouling can be reduced by using the appropriate shape of plate embossments or appropriate coatings [
3,
4,
5,
9,
10] or by choosing the optimal fluid velocity [
11,
12].
The optimisation of the cleaning process for PHEs at district heating substations could be aided by either the capability to determine the rate at which the layer of deposits builds up on the plates or by monitoring the thermal resistance of the heat exchangers. One of the methods proposed in the literature is to predict the formation of a deposit layer using numerical models of PHEs [
7,
8,
11]. All numerical models to predict the formation of deposits in PHEs, described in [
7,
8,
11], also considered their partial washout by the flowing fluid. In [
11], the authors have numerically analysed the formation of a layer of fouling in a plate heat exchanger using 3D Realisable κ-ε with non-equilibrium wall functions modelling and the von Karman analogy for three different apparatus (with different geometries) and taking into consideration the influence of geometry on fluid dynamics. The analyses were supported by experimental studies. The accuracy of the proposed solution is satisfactory, but its limitation is that the research can only be performed for turbulent flow. In addition, more complex models based on Computational Fluid Mechanics (CFD) tend to have long calculation times. Most articles investigating fouling in PHEs are related to the dairy industry, and articles [
7,
8] are examples of this. Both articles take a different approach to modelling the formation of the contamination layer from the previous one, based on the principle of conservation of mass and energy. In paper [
7], a 2D dynamic model of a PHE is proposed, allowing simulation with arbitrary cold and hot fluid flow configurations and with different fouling and cleaning kinetics. On the other hand, the paper [
8] discusses the problem of milk pasteurisation in a High-Temperature-Pasteurisation-Process plate pasteuriser. A model of the entire plant, including PHEs, was made. The study adopted a simple 1D model of the rate of fouling of the PHEs with protein with minimum calculation time and investigated the effect of the thickness of the resulting fouling layer on the milk pasteurisation temperature. Both approaches presented are not applicable to PHEs in district heating substations because of the different mechanisms of fouling formation due to the different types of working fluid.
Regardless of the increasing development of CFD models, the best and simplest method of evaluating contamination in PHEs is to determine their thermal resistance via experiments. An example of just such an approach is presented in [
13]. The authors of this paper aimed to determine the thermal resistance of a plate heat exchanger included in a district heating substation. Like the method developed and presented in this paper, the technique described by S.B. Genica et al. in [
13] compares the thermal resistance of a fouled and clean heat exchanger. However, in their study, based on the four plate heat exchanger models analysed, the thermal resistance of the fouling was determined as a function of shear stress. In addition, the overall heat transfer coefficient for a clean heat exchanger was determined based on its manufacturer’s design software. The need to use this software is a limitation of this method. A literature review indicates that there is still a need to develop a quick and simple method to optimise the PHEs cleaning process.
This paper proposes a method for determining the thermal resistance of fouled PHEs using experimental studies. It consists of determining the thermal resistance of a clean and a fouled heat exchanger. If the difference between these resistances exceeds a preset value, the heat exchanger will require cleaning. The thermal resistance of a fouled heat exchanger is determined by measurements from the transformed Peclet equation. At the same time, the thermal resistance of a clean heat exchanger is determined as the sum of the thermal resistances of heat transfer and heat conduction, where the thermal transfer resistances are counted from correlations on the Nusselt number of the cold and hot fluid. These, in turn, are calculated from the proposed modified iterative Wilson method. This paper presents the results of calculations using a new modification of the Wilson method on a constructed experimental test stand.
3. Determination of the Overall Heat Transfer Coefficient in the Heat Exchanger of a District Heating Substation
To determine the overall heat transfer coefficient
of the heat exchanger, the formula obtained from the Peclet equation can be used [
14]:
where
is the heat exchanger heat power in W,
is the correction factor that depends on the heat exchanger type (
= 1 for parallel-flow and counter-flow heat exchangers [
15]),
is overall heat transfer coefficient in W/(m
2K),
is the total heat transfer surface area in m
2, and
is the logarithmic mean temperature difference in K.
The heat output of the heat exchanger can be calculated from the energy balance equation for the hot fluid flowing through the apparatus, as , based on the known mass flow rate and its average specific heat , as well as the fluid temperatures at inlet and outlet .
Similarly, the heat flow rate in the heat exchanger can be calculated from the energy balance equation for the cold fluid, as , based on the known mass flow rate , average specific heat , and temperatures of the cold fluid at the inlet and outlet of the apparatus.
If there is no heat loss from the heat exchanger to the surroundings, then
. In practice, even though thermal insulation is used on the heat exchangers, heat losses occur and the heat flow rate value
is higher than the value
. In this case, the average value of the heat flow rate
that is transferred in the apparatus from hot to cold fluid can be calculated as [
16]
Based on the average heat flow rate
, the heat transfer surface area
and the mean logarithmic temperature difference
, the overall heat transfer coefficient
for the fouled heat exchangers can be obtained using the transformed Equation (5):
4. Determination of the Nusselt Number Correlation for Cold and Hot Fluid in a Clean Heat Exchanger
The method proposed in the article is experimental and modifies the popular Wilson graphical method. Thus, the laboratory studies to identify the cold fluid Nusselt number correlation are based on measuring cold and hot fluid parameters at a constant hot flow rate and for various cold fluid flow rate values. In turn, to determine the hot fluid’s Nusselt number correlation, laboratory testing in the heat exchanger must be carried out at a constant hot fluid flow rate and a variable flow rate of the hot fluid.
The developed method is based on iterative calculations. The single-iteration calculation algorithm presented in this chapter is for a single measurement point. Performing a series of iterations for several measurement points based on the described algorithm allows us to determine the Nusselt number correlation of the cold fluid .
At the beginning of the iteration, the Nusselt number
is calculated for the selected
measurement points using a modified Colburn correlation [
17]:
in which the constant
replaces the value of 0.023 and an additional constant
is inserted. In Equation (8),
and
denote the Reynolds and the Prandtl numbers of the cold fluid, respectively.
Based on the calculated Nusselt number
, the heat transfer coefficient at the surface washed by the cold fluid
is calculated using the following equation:
where
is the hydraulic diameter of the channel the cold fluid flows through in m and
is the cold fluid thermal conductivity coefficient in W/(mK).
Next, using the heat transfer coefficient on the cold side of the heat exchanger
, the heat transfer coefficient on the surface washed by the hot fluid
is determined using the transformed Equation (3) for the overall thermal resistance of the PHE:
Thereafter, based on the computed hot-side heat transfer coefficient
, we can compute the hot fluid Nusselt number
using the following equation:
where
is the hydraulic diameter of the channel the hot fluid flows through in m and
is the hot fluid thermal conductivity coefficient in W/(mK).
The calculated values of the Nusselt number
for
measurement points allow for the determination of the constants
and
via approximation with a linear function of the following form:
using the linear regression method.
In Equation (12), and are the Reynolds and the Prandtl numbers of hot fluid, respectively.
With the determined coefficients and (slope and vertical intercept of the function, respectively), the Nusselt number is recalculated using Equation (12).
In the next step of the algorithm, the hot-side heat transfer coefficient
is computed:
Based on the calculated heat transfer coefficient on the hot side of the heat exchanger
, the cold-side heat transfer coefficient
can be calculated as follows:
The Nusselt number of the cold fluid
is then determined as follows:
From the calculated Nusselt number and the Prandtl and Reynolds numbers of the cold fluid at the end of the iterative loop, the constants and are determined for measurement points by approximating the data with a linear function (8) using the linear regression method.
The iterations are continued as long as the following conditions are met:
and
where
j is the number of the iteration and
and
have small values (less than 0.001 for
and 0.01 for
)
The calculation algorithm for determining the constants and in the Nusselt number correlation for hot fluid can be written analogously to the algorithm for determining the constants and in the Nusselt number correlation for cold fluid presented above.
In the single-iteration loop, an apostrophe “ ” has been used in some symbols, i.e., in the heat transfer coefficient and the Nusselt number for cold and hot fluid. This was performed as such because, in the single-iteration loop, these values appear and are calculated twice. The rule of thumb was adopted, with the symbol “ ” denoting the heat transfer coefficients calculated from the formula for the overall heat transfer resistance between the cold and hot fluid in the heat exchanger (Equations (10) and (14)) and the Nusselt numbers calculated from them (Equations (11) and (15)).
Plate heat exchangers in district heating substations may differ in design, i.e., in the number of plates or the type of corrugations. The proposed computational algorithm will allow the determination of constants in the correlations on the Nusselt numbers for cold and hot fluid flowing through a specific clean apparatus model. This, in turn, will make it possible to determine its total thermal resistance for any operation conditions.
5. Description of Laboratory Stands for the Verification of the Modified Wilson Method
The research discussed in this article will be conducted at two laboratory stands. One is a district heating substation belonging to the Municipal Heat Engineering Company (MPEC) in Krakow, which has existed for many years in one of the Cracow University of Technology buildings (
Figure 1). The other is a new laboratory stand, a heat substation located in the laboratory of the Department of Energy (DE) of the Cracow University of Technology (
Figure 2 and
Figure 3a).
The heating substation in the DE laboratory is a parallel substation supplying two heating systems: a domestic hot water system (with PHE 1 in
Figure 2) and a central heating system (with PHE 2 in
Figure 2). The substation is supplied with hot water from an individual heat source, an electric boiler of the company ELTERM of 15 kW. The hot water reservoir used during the tests consists of two storage tanks with a capacity of 800 litres. A Danfoss heat exchanger model XB12M-1-16 (composed of 16 plates) is installed in the central heating system. In comparison, a Hexonic heat exchanger model LC110-30-2 (consisting of 30 plates) is installed in the domestic hot water system (
Figure 3a). A drawing of a single heat exchanger plate with dimensions is shown in
Figure 3b, while
Figure 3c shows the type of embossing on the plates (chevron). Both heat exchangers are used as counter-flows.
The flow rate of cold and hot fluid in the system is measured using axial turbine flow meters ranging from 4 to 160 L/min (measuring accuracy: ±3% of measured value). Water temperatures in the installation are measured using type K (NiCr-Ni) sheathed thermocouple sensors with an outer diameter of 1.5 mm class 1.
Figure 2 shows the location of sensors for measuring the temperature and flow rate. The location of sensors in the hot water system with PHE 1 are as follows: Th1,1 and Th2,1—hot water temperature measurement at the inlet and outlet, Tc1,1 and Tc2,1—cold water temperature measurement at the inlet and outlet, and mc1 and mh1—cold and hot water flow rate measurement. The location of sensors in the central heating system with PHE 2 are as follows: Th1,2 and Th2,2—measure the hot water temperature at the inlet and outlet, Tc1,2 and Tc2,2—measure the cold water temperature at the inlet and outlet, and mc2 and mh2—measure the cold and hot water flow rate. Compact pressure sensors with piezo-resistive measuring cells with temperature compensation and a pressure membrane made of special steel, whose basic error is ±0.1% of the final value, are used to measure the pressure in the installation. In contrast, the pressure drop in the heat exchangers is measured using differential pressure transmitters whose error for the basic measuring range is less than 0.05%. All sensors used for the measurements are from Ahlborn. All sensor measurement data are recorded by Ahlborn’s ALMEMO logger. The ALMEMO system enables real-time visualisation of the data collected from the sensors, data logging and the generation of reports in tables and graphs. Flow rate control in both systems is achieved via control valves. The research described in detail in this article will be carried out at the second stand.
6. Computation of the Cold Fluid Heat Transfer Coefficient in the Clean Heat Exchanger Using Measurement Data
The proposed algorithm, which is a modification of the Wilson method, was verified by carrying out test computations. They used measured data from the heat exchanger in the domestic hot water system of the laboratory stand described in the previous section.
Measurements carried out during the laboratory tests included measurements of water temperature at the inlet and outlet of the heat exchanger, as well as flow rates of cold and hot water. The data sampling time during the measurements was 1 s. Cold and hot water physical properties were determined for the average temperature values at the inlet and outlet of the heat exchanger, separately for each fluid.
The domestic hot water system under study uses a Hexonic LC110-30-2 PHE, and its heat transfer area . The correction factor , because the analysed heat exchanger is a counter-flow one.
The measurements were carried out so that the flow rate of cold water was constant (about 40 L/min.) while the flow rate of hot water was changed (in the range from about 5 to 40 L/min). From the archived database,
measurement points were selected. At the beginning of the test calculations, the mean heat flow rates
and overall heat transfer coefficients
of the PHE under test were determined for a selected number of measurement points (
). The results of these calculations are shown in
Table 1.
Table 2 lists the average values of the temperatures
and
and the water properties, i.e., densities
and
, kinematic viscosities
and
, Prandtl numbers
and
and heat transfer coefficients
and
, for hot and cold water, respectively.
In the next step, the cold water Nusselt numbers
were calculated according to the modified Colburn correlation expressed by Equation (8). The initial values of the constants
and
in the first iterative loop for calculating the Nusselt number values
were taken as
and
. The water’s physical properties were computed for the average temperature using [
18]. The calculated values of the Nusselt number for the selected measurement points
are shown in
Table 3. In
Table 3 et seq.,
denotes the number of measurement points in the sequence.
Then, the cold-side heat transfer coefficients of the heat exchanger were calculated using Equation (9), and later, the hot-side heat transfer coefficients were calculated using Equation (10). In order to calculate , it is necessary to know the thickness of the plates in the heat exchanger, which is , and the thermal conductivity of the material from which they are made (stainless steel), which is . In the next step, based on Equation (11), the values of the Nusselt number for hot water were determined, which were then approximated by the linear function (12), which allowed us to determine the constants and , the values of which are −0.23977 and 40.00581, respectively.
The results of the subsequent calculations listed above needed to calculate the constants
and
are shown in
Table 4.
Table 4 also contains the values of the Reynolds
and Prandtl
numbers of hot water necessary to perform the approximation with a linear function.
In
Table 4, for the first measurement point, the heat transfer coefficient
and the Nusselt number
have negative values. This only occurs in the first iteration and is due to the fact that when determining
from Equation (10), the sum of the heat transfer resistance on the cold side and heat conduction was larger than the total heat transfer resistance of the heat exchanger. Therefore,
in this case has a negative value. The Nusselt number
was determined from the
of Equation (11), so its value is also negative. The problem was eliminated in subsequent iterations, due to the increasingly better fit of the constants in the Nusselt number correlations.
After determining the constants
and
, the Nusselt numbers
from Equation (12) were again determined, and from them, the heat transfer coefficients
from Equation (13) were calculated, followed by the heat transfer coefficients
from Equation (14). The Nusselt number values for cold water
were then computed from Equation (15). The results of the calculations are presented in
Table 5.
The first iteration ends with the approximation of the calculated values of the Nusselt numbers
by Equation (8) using the linear regression method, resulting in the constants
and
, which are 0.01988 and 2.55663, respectively. The coefficient of determination of the approximation using the linear function (8) performed is
. The approximation using the linear function (8) is also shown graphically in
Figure 4. The computed constants
and
are the initial values for the next iteration loop.
The computational algorithm for determining the Nusselt number correlation for a cold fluid is shown as a block diagram in
Figure 5.
Table 6 shows the calculated Nusselt numbers
and heat transfer coefficients
for cold water in five successive iteration steps.
The computational results of the successive iterations presented in
Table 6 show that the proposed iterative method converges. The values of the constants
and
determined in subsequent iterations and the coefficients of determination
are summarised in
Table 7 and shown graphically in
Figure 6a,b.
The modified Wilson method described in this article enables the determination of constants in the correlation of the Nusselt number for cold fluid flowing through a heat exchanger. Both conditions (16) and (17) were satisfied after 17 iterations. The final form of the correlation for the Nusselt number for cold water, when the hot water flow rate is about 40 L/min, is as follows:
However, it can also be used to determine the constants in the correlation of the Nusselt number of a hot fluid, but this will require subsequent measurements in which the hot fluid flow is constant and the cold fluid flow is variable.