1. Introduction
Changing the direction of thermal fluid flow by 180° significantly expands the design options for flow devices. One of the most common solutions used for this purpose is the use of a u-shaped elbow (also known as return bend or u-tube). Such a solution makes it possible to reduce the size of apparatuses while preventing the problems of elongation compensation associated with thermal expansion of tubular elements, which, especially for high-temperature thermal fluids, is common. One apparatus where the u-tube elbow design solution has been widely adopted is the compact shell-and-tube heat exchanger. Each tube bundle of such an exchanger is terminated on one side by a u-tube element that returns the thermal fluid flowing in the tube toward a head mounted on the shell and equipped with both inlet and outlet ports. Such a solution, due to the very favorable relationship of achieved parameters to external dimensions, is often used in industrial equipment and installations. They are encountered, among others, as devices in compressor systems, combustion engines and chemical reactors [
1,
2]. In addition, u-tubes themselves are often used for heat exchange purposes. This is the case, for example, in solar collectors [
3,
4], ground heat exchangers [
5,
6,
7], and surface cooling systems [
3,
8,
9].
Extensive use in engineering would not be possible without thorough research into the nature of flow in such elements. The complex and dynamic flow phenomena occurring in u-tubes directly affect the distribution and velocity profile of the flowing thermal medium and thus the performance of the implemented processes. It is of significance that the nature of flow in the u-tube intensifies the erosion and corrosion processes on the inner surfaces of the pipe walls [
10,
11,
12]. A comprehensive analysis of these issues is discussed in the paper [
13]. This is a significant operational complication of u-tube-based flow systems. It should also be noted, in most available studies, that flow phenomena occurring in tube bends are combined with the flow velocity of the medium and the bending radius of the u-tube. For this reason, in the case of compact shell-and-tube heat exchangers where, most often within a single unit, we have to deal with different geometries of the u-tube tube bundle, these construction–design issues become sensitive [
4]. One of the first areas of deep study of u-tube-equipped apparatus was strength analyses related to stresses in u-tube support baffles [
14] or with heat transfer modeling [
15]. The Finite Element Method (FEM) has been widely used for this purpose for years. On the other hand, in the area of studying flow phenomena occurring in curved geometries, the dominant method is Computational Fluid Dynamics (CFD) and using genetic algorithms [
16,
17], as well as experimental techniques, including Particle Image Velocimetry (PIV) [
18]. It should also be noted that most of the available studies are in the recommended range of fluid flow velocities and bending radii are clearly larger than the diameter of the tube being tested. For example, the paper [
19] showed that the u-tube element has little effect on the heat transfer conditions of the straight section of pipe in front of it, while it significantly shapes the flow conditions behind the u-tube element. For the conditions studied, this effect terminated no sooner than at a distance of 90d behind the u-tube. In turn, in the paper [
19] it was noted that with an increase in turbulence generated by the bending radius, diameter of the u-tube or flow rate, the enhancement of heat transfer efficiency does not exceed the value of 25% in relation to the straight pipe, and is quite limited in range and quickly disappear in the continued flow through the pipe. Under laminar flow conditions, a different correlation was noted. Studies [
20] have shown that in such a regime the flow of the thermal fluid is disturbed to a much greater extent, and the heat transfer coefficient under such conditions can increase by up to 3–5 times.
The literature also contains a number of papers on two-phase gas–liquid flow in u-tubes [
21,
22,
23]. However, the implementation of the results obtained in the area of research undertaken in this paper was limited.
One of the few works where u-tubes with small bending radius were studied is the work of [
24]. It focused, among other things, on the analysis of Dean and Reynolds numbers in the knee region and in straight sections. It was found that, despite the fact that in the elbow area the Dean number has a greater influence on heat transfer than the Reynolds number, the Reynolds number should be based on the design calculations due to the larger area of straight sections. The second conclusion of this work, which is important from the point of view of the presented research, is that turbulence intensity is strongly related to improvements in heat transfer calculations in such systems. However, work such as the above and others in the field of u-tube hydrodynamics do not address the range of geometric and flow parameters studied in this paper that are typical of compact heat exchanger tube bundles.
For such flow systems, no information has been found regarding the maldistribution of fluid in the straight region of the pipe behind the elbow. This is a barrier for heat exchanger manufacturers to undertake design work on new u-tube bundle layouts. The strong dependence of the maldistribution of flow in the u-tube on the flow conditions of the thermal fluid makes it difficult in the application tasks to carry out the selection of the geometry and flow parameters. Tests and measurements taken under conditions as close to real as possible are helpful in this field. For this reason, it was decided to pursue the concept of testing the geometry of u-tubes used in a real compact heat exchanger operating at velocities slightly above the values accepted as economical for a heat exchanger as an interesting alternative to increase the performance of the apparatus within the still-acceptable increases in pressure drop. Thus, the purpose of the ongoing research is to parameterize the maldistribution of flow under conditions appropriate to real tube bundles of compact heat exchangers, but importantly at increased flow velocities. Based on the available literature, it was found that, in addition to the CFD method commonly used for this purpose, the study will be extended to non-invasive tests using an optical method. The results obtained in this way will extend the description of velocity distributions and allow analysis of the intensity of flow turbulence in u-tube systems.
2. Materials and Methods
The research methodology consisted of two main parts: simulation studies using the CFD method and experimental studies using the PIV technique. Both parts of the study were conducted under analogous geometric and flow conditions. The tests adopted pipe models with internal diameter d = 0.007 m, length l = 0.205 m, and bending radii r
g of: 0.039 m; 0.023 m; 0.009 m, resulting in d/r
g ratios of 0.18, 0.3, and 0.78, respectively. The experimental models were made of transparent poly(methyl methacrylate) material PMMA. The liquid supply was provided through the upper port of the u-tube and the liquid outlet occurred in its lower port. The schematic diagram of the test rig is presented in
Figure 1. Water with the addition of inert fluorescent markers (FPP-RhB-10), based on the melamine resin necessary for the PIV technique, was used as the fluid in the tests. The tests were conducted for four liquid Q rates of: 0.2 m
3/h, 0.3 m
3/h, 0.4 m
3/h, 0.5 m
3/h, and corresponding to superficial liquid velocities v
p equal to 1.44 m/s, 2.17 m/s, 2.89 m/s, 3.61 m/s, respectively. The first two values corresponded to the recommended and rate limit of liquid flow in the assumed geometry, and the next two correspond to the increased flow velocities. The total number of simulation and experiment variants was 24.
Table 1, in turn, presents a summary of the main geometry and flow parameters during the conducted tests.
2.1. CFD Metodology
In simulation studies using CFD, the evaluation of the distribution of velocity and turbulence intensity of the flowing medium was carried out on the basis of the finite volume method based on the code of the Ansys Fluent software version 19.2. Due to the fact that there is turbulent flow in the analyzed computational domain, there was a need for its faithful representation. Several methods for simulating turbulent phenomena are present in the literature. One can find a considerable number of references to methods: Direct Numerical Simulations (DNS), Reynolds-Averaged Navier–Stokes Equations (RANS), and Large Eddy Simulation (LES) [
25]. However, the most widely used is the RANS method. It comes under several variations, based on the following models: the Reynolds Stress Model (RSM), the Boussinesq models (k-ɛ and k-ω), the Prandtl model, and the Kolmogorov model. The k-ɛ model was chosen for the study. The use of this model ensures a high correspondence between CFD results and experimental studies with the minimum possible load on computational units. This model has been successfully used in CFD studies for flow in heat exchangers [
26,
27,
28,
29,
30,
31].
The main equations of this method in the computational domain can be rearranged in the following figures [
32]:
Energy equation:
where:
u—averaged velocity of the fluid [m/s];
p—pressure;
ρ—density of the fluid;
T—temperature;
µ—kinematic viscosity of the fluid [m
2/s];
Cp—specific heat capacity [J/kgK];
λ—thermal conductivity [W/mK].
On the other hand, the k-ɛ model itself can be presented as follows [
33]:
Kinetic energy of turbulence:
where:
where:
k is turbulent kinetic energy [m
2/s
2],
ɛ is turbulent dissipation rate [m
2/s
2],
Gk is producing term of turbulent kinetic energy generated by mean velocity gradient, C
1ɛ and C
2ɛ are constants, σ
ɛ and σ
k are Prandtl numbers corresponding to turbulent kinetic energy and turbulent dissipation rate, and
µt [Pas] is expressed as:
where:
Cµ = 0.09 [-], C
1ɛ = 1.44 [-], C
2ɛ = 1.92 [-], σ
k = 1.0 [-], σ
ɛ = 1.3 [-], and
Gk [-] is defined as:
The study was conducted in a transient state (time step was 0.00001). The pressure-velocity coupling was performed with the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm with the PRESTO (PREssure STaggered Option) pressure discretization scheme. Second-order upwind interpolation was used to determine representative samples of component values on the surface of control volumes and a standard wall function. Variable values under relaxation factors were determined. They were, respectively, 0.6, 0.5, 0.8, and 0.8 for the pressure, momentum, turbulent kinetic energy, and turbulent energy dissipation, respectively. The following convergence criteria were adopted: 1 × 10−7 for the continuity equations, 1 × 10−5 for others. Boundary conditions for the inlet: ‘velocity inlet’ (turbulence intensity: 5%, hydraulic diameter: 0.007); for the outlet: ‘pressure outlet’.
The study was conducted for the 3D computational domain. Three geometries were generated, corresponding to the bending radii of the u-tube shown in
Table 1. Due to the fact that a key element of CFD studies is the discretization of the computational domain, an analysis of the sensitivity of the mesh density to the results obtained was carried out. This study was conducted for a u-tube with a bending radius of r
g = 0.023 m (
Figure 2; the figure additionally includes the location of the planes for which the analysis was conducted in
Section 3). A variable-density (hexagonal) grid was used. Four grid densities ranging from about 4 × 10
5 to about 1.2 × 10
6 elements were generated. Sensitivity analysis was carried out in terms of average in-plane velocities in the axis of symmetry of the elbow (for Q = 0.5 m
3/h). The values of average velocities were, respectively: grid 1 (401,289 elements): 3.772 m/s; grid 2 (521,487 elements): 3.711 m/s; grid 3 (606,732 elements): 3.629 m/s; grid 4 (1,119,562 elements): 3.621 m/s. The difference between the minimum and maximum values was about 4%. It was found that an increase in grid density did not have a large effect on the obtained results. In addition,
Figure 3 presents the velocity profiles at the axis of symmetry of the elbow (the exact location is marked in
Figure 3 with a yellow line) for the four grid densities. It was found that the profiles for grids 3 and 4 follow a similar trend, and there are few differences. Therefore, the grid density corresponding to grid 3 parameters was used for further studies (
Figure 4). This reduced the time consumption of the CFD testing process. For the remaining u-tube with different bending radii, discretization parameters corresponding to grid 3 density were assigned.
The next task was to validate the adopted assumptions of the numerical studies and the process of discretization of the computational domain. Validation of the CFD results was carried out with reference to the obtained values of pressure drop during the experimental tests and velocity profiles in the axis of symmetry of the elbow, with reference to the PIV results. As in the case of sensitivity analysis of the computational grid, the validation process was carried out for a u-tube with a bending radius r
g = 0.023 m and a flow rate Q = 0.5 m
3/h. The obtained value of pressure drop by the CFD method was 12,550 Pa. The calculation error was about 6% (the value obtained by the experimental method was 11,840 Pa). In turn,
Figure 5 presents the velocity profiles obtained by the CFD method and the PIV method. It can be noted that the velocity profile obtained by the CFD method is close to the reference values obtained by the PIV method. In addition, the average velocity values in this profile for both test methods were compared. The average velocity obtained by the CFD method was 2.64 m/s. In contrast, the value for PIV testing was 2.98 m/s. The difference between the two methods was less than 12%. The presented data confirm the correct parameterization of the assumptions of numerical tests.
The observed differences between the CFD and PIV methods may have been due to some limitations of both methods. Numerical models have some restrictions or simplifications. As a consequence, the obtained results are characterized by greater regularity and stability. In the case of PIV studies, on the other hand, the accuracy of the results may be affected by the optical inadequacy of the material of which the model is made (this is especially valid for the wall region), factors related to the apparatus used (the type of camera, its focus, reflection of the laser light), and the human factor. The uncertainty of PIV measurement is interpreted more extensively in
Section 2.2.
2.2. PIV Metodology
Experimental tests following the numerical calculations were carried out using the PIV method for analogous geometric and flow assumptions as in the numerical studies (
Table 1). The u-tube models were cored from transparent PMMA material. In order to avoid barrel distortion generated by curved surfaces [
34], the u-tube models had a rectangular cross-section, and the execution in them by the cavity method of an internal channel with a diameter of d = 0.007 m made it possible to realize the flow as in a classical u-tube. A similar solution was presented in the work [
35]. Liquid supply was carried out through the upper port of the u-shaped element and liquid outlet occurred in its lower port. The stubs were connected to the pumping system and the tank so that the movement of the liquid took place in a closed circuit. An ENKO MPP-05C electromagnetic flow meter was responsible for measuring the fluid flow. In addition, a Peltron PXWD-0.2 differential pressure transducer was used to measure the pressure drop between the inlet and outlet ports of the test unit. The temperature of the liquid during the experimental study was measured and varied between 20.1 and 20.4 °C (inlet of the tank). Thus, it was assumed that there was no temperature gradient that could affect changes in the nature of the fluid flow. A constant temperature value in the laboratory and a fluid distribution system with a large volume (190 dm
3) in relation to the volume of circulating fluid were certainly responsible for this situation. The optical path of the test rig was formed by a Dantec Dynamics FlowSense EO 4M CCD camera positioned perpendicular to the test object and the laser light sheet generated by the Dantec Dynamics DualPower TR laser (optical knife technique). The camera was equipped with an Omega Optical 550LP narrow-band filter to facilitate imaging of light generated by fluorescent markers. A Berkeley Nucleonics Corp model 575-8 pulse generator was responsible for synchronizing the camera and laser. Control of the equipment and acquisition, data processing, and performance of the necessary PIV calculations were carried out in the Dantec Dynamics Dynamic Studio environment.
Each u-tube element of the tube bundle was studied separately, in the same geometric position, and for all assumed flow parameters. Images were recorded at 2048 × 2048 pixel resolution, in double frame mode at 10 Hz. The time between laser pulses was selected using the PIV Setup Assist tool provided by the manufacturer of the PIV system and was set at 100 μs to meet the requirements of the measurement method [
36]. Each measurement series included the registration of 100 double frame images. The study thus included 12 measurement series. The first stage of the PIV analysis was image processing (selection of the measurement area, balancing the brightness of the images, determining the time-averaged image, and extracting moving markers using arithmetic calculations). The obtained image processing results made it possible to determine the vector velocity field based on the Adaptive PIV algorithm. The obtained vector velocity fields were used to determine scalar maps of velocity values V defined as
and scalar maps of turbulent flow intensity factor
Ti defined as
similarly as in the work [
37]. Details of the PIV system are summarized in the
Table 2.
An important part of any measurement method is the determination of measurement accuracy. The PIV method is complex in this aspect, however, based on the assumptions presented in the work [
36], it can be stated that the most significant contribution to this aspect is the estimation of the seeding particles displacement uncertainty σ, as this parameter directly affects the value of the found velocity vectors. To better illustrate the uncertainty σ, it is directly converted into a unit of velocity. The uncertainty distributions obtained during the study are shown in
Figure 6. They were determined based on the Peak Height Ratio as proposed in the work [
38], and, taking into account the cumulative Rayleigh distribution function, it can be found that 95% confidence was obtained that the measured velocity differs from the true value by, at most, σ. The uncertainty σ increases as the flow velocity increases and reaches a maximum value of 0.269 m/s for the elbow with the largest bending radius, which is due to the largest velocity values measured in this geometry.
3. Results
Based on the results of the flow visualization seen in
Figure 7, the separation of two areas with different average fluid velocities was identified across the measurement section of the u-tube, both in the bend section and in the straight section (before and after the elbow). Depending on the particular section of the u-tube (a bend or a straight section), the spatial pattern of liquid velocity distribution in these sections took on a different characteristic form. Studying the maldistribution in the straight section after the bend of the u-tube, it was shown that for the geometry of the u-tube model with r
g = 0.039 m (d/r
g = 0.18), the degree of flow maldistribution practically assumed a constant value, thus the maldistribution of velocity affected the entire analyzed straight section. The difference in velocity under such conditions, between the inner (upper) and outer (lower) streams, reached up to 25%. Higher velocity values were recorded in the lower part of the straight section of the u-tube. The observation of such a significant influence of the elbow under these flow conditions confirms the observations of the work [
20], where for stabilized, undisturbed flows the influence of the elbow was seen at significant distances of the straight section behind the elbow.
The presented study also showed that the process of stabilizing the flow in the straight section of the u-tube increased in intensity as the bending radius decreased. This relationship is a consequence of the introduction of additional turbulence into the flow of the thermal fluid as a result of the decrease in the bending radius. The shortest distance at which the return to steady-state flow conditions occurred was observed for geometry rg = 0.009 m, with velocity vp = 1.44 m/s. This distance was 17d, counting from the end of the bend. Such a value should be considered significant. Therefore, the simplified calculation method based on the averaged flow velocity of thermal fluid that contacts the tube surface uniformly in terms of velocity should be verified. Thus, it is suggested to determine the surface fractions of low- and high-velocity flow streams and perform thermal calculations separately.
When analyzing the local nature of the flow, the focus was on evaluating the velocity profile in the axis of symmetry of the bend section of the u-tube. There was a characteristic profile distribution in all analyzed variants. During the analysis, a strong correlation of the velocity profile distribution with the r
g value was detected. The localization of the areas with the highest and lowest fluid velocity for extreme values of r
g was reversed (
Figure 8). For the largest bending radius, the maximum velocity was recorded at a distance of 0.5 mm from the outer wall, while for the smallest bending radius the maximum velocity was recorded in the area closer to the opposite wall of the bend. Counting from the outer wall of the bend, this distance was 5.5 mm. The distribution of fluid flow velocity depending on r
g is characterized by the spatial inversion of fluid streams, which is close to a symmetrical distribution. This fact should be noted, especially in the aspect of evaluating the erosive action of the flowing fluid. Thus, the location of areas threatened by an increase in the intensity of erosive processes is directly related to the position of the u-tube in the tube bundle.
In the literature, one can find papers indicating a significant change in the nature of the flow depending on the imaging plane. Based on the observations of Mariotti et al. [
39,
40], visualization of the flow in several planes was carried out (see
Figure 2). The CFD method was used. For the results obtained along the u-tube (
Figure 9 and
Figure 10), clear fluid maldistribution structures were observed in cross-sections. The maldistribution phenomena are permanent and hold with the course of the u-tube (developing and disappearing depending on the flow velocity and bending radius). This confirms previous observations indicating that the symmetrical distribution of velocity in a straight inlet section is disturbed by the presence of a bend in every case analyzed. It should be considered as a point that initiates the movement of a stream with increased flow velocity to the outside of the elbow in the case of bending radii r
g of 0.039 m and 0.023 m and to the inside of the elbow in the case of bending radius r
g = 0.009 m.
In investigating the flow of thermal fluid inside a u-tube bundle, in addition to knowledge of the velocity distribution of the medium, information about the intensity of its turbulence is also important. A parameter that well-describes the quality of turbulence and at the same time can be determined by direct measurement of the component velocities is the turbulent flow intensity factor T
i. Analyzing the obtained test results, the specific effect of the bending radius on this parameter was noted (
Figure 11).
In the case of r
g = 0.039 m and r
g = 0.009 m, relatively small areas of increased T
i parameter were identified. In contrast, in the case of r
g = 0.023 m, the turbulence described by T
i parameter took on increased values over a much larger area. It was found that hydrodynamic conditions inhibited the formation of fluctuations in the pattern of the thermal fluid flow and were responsible for the reduced values of the T
i parameter. In the case of r
g = 0.009 m, such a factor was the abrupt change in the direction of the flow. In the case of r
g = 0.039 m, the stabilized separation of the fluid stream at the bend was that factor. Consequently, fluctuations in the velocity components u and v, which define the parameter T
i, were limited. In contrast, in the case of u-tube with r
g = 0.023 m, both the factors discussed above weakened. Thus, conditions were created, making the fluctuation of the u- and v-component velocities possible over an increasing area, which directly yielded the effect of increasing the value of the T
i parameter. This is also confirmed by the irregular pattern of the edge of the area where velocity V is increased, visible in
Figure 6, for u-tube with an average bending radius r
g.
In order to better understand the mechanism of the formation of fluid flow maldistribution in the tested u-tubes, streamlines were determined. It was noted that the paths of fluid movement are strictly determined by the geometry of the pipe bending. In each case studied, similar flow patterns were obtained. The fluid stream characterized by a higher velocity moving along the elbow comes closer to its outer wall (fluid stream drift). As a result, a maldistribution of streams with different velocities was formed at the outlet of the elbow. Visualization of the approach of the higher velocity fluid stream to the wall makes it possible to identify the areas most sensitive to the erosive action of the fluid. The identified flow asymmetry increases as the bending radius decreases. At the same time, for u-tubes with rg of 0.039 m and 0.023 m, it is significantly less identifiable. The strongest drift of the fluid stream toward the outer wall of the elbow occurs at a bending radius of rg = 0.009 m. On the other hand, analyzing the effect of the change in fluid velocity on the asymmetry of the streamline, it was noted that the asymmetry increases as the fluid flow rate increases. This confirms the strong dependence of the maldistribution of the fluid in the u-tube on these two parameters. Protection against excessive erosion wear of u-tubes is the main physical aspect from the study. This is the most serious operational problem because replacement of the inner elbow of a u-tube bundle is impossible without removing the outer bundles, which raises the cost of such an operation.
It should also be noted that no turbulence or secondary flow patterns were observed in the studied longitudinal section of the u-tube.
Figure 12 presents a case of the most evident of the observed streamline asymmetries occurring with fluid flow through an elbow of r
g 0.009 m. In other cases, the differences in streamlines distributions are hardly visible with the range of parameters adopted in the tests.
Introducing the relationship between d/r
g and the T
i parameter into the description of the hydrodynamics of thermal fluid in u-tube systems, the course of changes in the maldistribution of thermal fluid in the studied u-tube geometries was determined (
Figure 13). By performing the above analysis as a function of the fluid flow velocity, it was shown that the turbulent flow intensity factor T
i for the recommended flow velocities is distributed differently compared to increased velocities. A significant decrease in T
i gain with increasing fluid velocity may result in a smaller heat transfer process efficiency increment than expected from the increased flow velocity. The limiting value of the geometric parameter d/r
g, above which a change in curvature no longer causes significant increases in flow turbulence, was found to be 0.3. Above this value, the increase in the T
i coefficient does not exceed 10.2% of its maximum value (the case of v
p = 2.29 m/s), and it should be remembered that the flow of the fluid under these conditions causes an increase in pressure drop of 125 Pa.
The findings from such a parameterization of the maldistribution of thermal fluid flow in u-tubes support design and construction decisions. Considering the pressure drop increment and the limit value of the d/rg parameter for the increase in flow turbulence, it is inappropriate to consider deliberately increasing the thermal fluid flow velocity in apparatuses of this type. In addition, due to the high surface fraction of areas with an inhomogeneous pattern both in terms of fluid velocity and spatial distribution, it is necessary to consider the possibility of varying the calculation of the heat transfer process along the length of the u-tube in the tube bundle of the heat exchanger. Taking into account the surface share in the velocity domain, at the stage of the thermal and flow calculations, one can increase the accuracy of the algorithms used for the selection of compact heat exchanger geometries. These issues appear to be particularly relevant for the case of external tubes in a bundle with large bending radii, where the rate of flow stabilization is low.
As written in the introduction, u-tubes are a common part of industrial flow infrastructure. Therefore, in any application other than heat exchangers, when the geometry and flow parameters are similar to those analyzed in the manuscript, similar fluid maldistribution should be expected. The typical way to proceed should be to determine an area of the u-tube that plays an important role in the process being implemented. If the area of interest of the analyzed process coincides with the area where maldistributions were identified in this study, the effect of stream partitioning due to flow velocity should also be considered in the procedural calculations. It is clear that such an approach is justified only for processes driven by velocity or its gradient.