Numerical Study of Flow Maldistribution in Multi-Plate Heat Exchangers Based on Robust 2D Model

: Plate heat exchangers (PHE) are characterized by high heat transfer efﬁciency and compactness. An exploitation problem of the PHE is related to ﬂow maldistribution, which can make part of the PHE idle, resulting in overheating and damage. Making geometrical modiﬁcations to the PHE can help reduce ﬂow maldistribution. Modiﬁcations should be kept to a minimum, so as not to complicate the production process. There is a large number of possible geometrical modiﬁcations, which simply considers additional obstacles or stream dividers. To test all of them would be impractical and would also take a prohibitively long amount of time to obtain experimental measurements. A typical PHE is characterized by a complex system of channels. Making numerical calculations of its 3D model can be prohibitively time and resource-consuming. The present work introduces a physically consistent methodology of the transformation of a real 3D geometry to its 2D representation. Its main novelty is to assure the same pressure drop balance remains between the 3D and 2D geometries. This is achieved by a preservation of the same cumulative pressure losses in both geometries. The proposed innovative approach levels the pressure balance difference by adding properly designed local geometrical modiﬁcations. The developed methodology allowed a wide range of parameter space and various geometrical modiﬁcations to be investigated, and revealed geometrical optimizations leading to the improved performance of the PHE. To minimize the inﬂuence of other factors, an incompressible and single-phase ﬂow was studied.


Introduction
The plate heat exchanger (PHE) has a very efficient design but when compared to other design solutions, can demonstrate an important exploitation problem due to uneven flow distribution in the channels in-between its plates. The experimental and numerical findings of numerous authors confirm the existence of a flow maldistribution problem in practically every type of PHE design. The work of Mueller et al. provides a comprehensive review of different possible causes that may lead to flow maldistribution [1]. This may be provoked by mechanically driven factors such as fouling, fabrication tolerances, a bypass or poor header performance, or can also be related to the thermodynamic character of the ongoing processes (two-phase instabilities, heat transfer induced by changes of viscosity or density) [2].
The flow maldistribution problem is reported in non-phase-change applications [3][4][5][6] as well as in PHE condensers [7][8][9] and evaporators [10][11][12][13]. The uneven load of the PHE channels is especially undesirable in the case of evaporators. It can create dead zones in the farthest channels on the refrigerant side of the PHE and the vapor may be trapped in certain parts of the device, causing a local dry out or overheating. Hence, the overall heat transfer efficiency of the PHE may decrease significantly.
Many solutions to decrease the magnitude of the flow maldistribution have been proposed in different subcategories of the PHE, and they are mainly focused on inlet header design modifications. In the case of multi-stream plate fin exchangers, solutions include additional inlet baffles [14,15], rearranged channel configurations [16] or the introduction of internal two-phase distributors to the inside of the exchanger [17]. The work [11] shows the possible improvements that may be made to the performance of a two-stream PHE by an inlet distributor modification.
The flow maldistribution in the two-stream PHE can be potentially improved by simply rearranging the inlet and outlet locations. When the streams enter the devices from opposite sides of the PHE (Z-type), the mass flow across the flow channels are shown to be more evenly distributed compared to that of a (U-type) design [18]. However, the Z-type solution is not commonly preferred in HVAC installations for practical reasons, e.g., to reduce occupied space, and for maintenance and conservation.
The current studies are built on the assumption that flow maldistribution can be minimized with geometrical optimizations to the PHE. The easiest modification includes the introduction of additional obstacles or separators to the header pipe of the PHE. In practice, there is a large number of possible modifications. To examine all of them experimentally is impractical and would be greatly time consuming. Typically, cases like this would benefit from a virtual experiment. The main challenge in the numerical modeling of a real PHE is related to its very complicated geometry and the requirement of a high-resolution mesh to properly capture the flow as it relates to physics and heat transfer. Consequently, exploring a wide range of flow parameters and geometrical variations can be even more time consuming, or even impossible, for averaged-sized supercomputers [19]. This problem can be solved by an appropriately designed calculation methodology based on physically consistent simplifications that lead to a fast and reliable numerical model.
In our previous work [19], the properly designed 3D numerical mesh of a typical PHE was shown to consist of hundreds of millions of computational nodes. This is prohibitively too much for engineering applications, especially for flows with a phase change. Moreover, it makes the investigation of a wide range of parameter space and the geometrical modifications of a heat exchanger impossible. The current available research, which uses 3D numerical calculations, is limited to a single pair of plates [20] or focused on a precise detail of a 3D flow feature [21]. Opposite to this, if full geometry is being considered, 2D numerical calculations are still exploited [22].
The main objective of the present study is to create a robust methodology to transform the 3D geometry of a real PHE to its 2D representation. Consequently, fast and physically consistent calculations of a flow in the PHE are also made possible. In the above-mentioned author's work [19], a reduction of the 3D geometry to its 2D representation has been successfully carried out and confirmed through comparison with experimental results. Nevertheless, the work only considered a single pair of plates (one flow channel). In the current studies, a multi-plate PHE is considered. This makes it possible to investigate the flow maldistribution phenomenon and different possible ways of reducing it.
The previously proposed methodology of the 2D reduction assured a preservation of geometrical and flow features, such as a channel Reynolds number and a mass flow. In the case of a multi-plate PHE, this is not enough because of significantly important differences in the pressure drop balance between the original 3D geometry and its 2D representation. The main novelty of the presently proposed numerical approach is in the development of a methodology assuring the same pressure drop balance between the 3D and 2D geometries. This is achieved by preserving the same cumulative pressure losses in both geometries. The author's original idea was to level this difference in pressure balance by adding a properly designed local geometrical modification. This required a large number of careful numerical calculations to establish a trustworthy relation of this local pressure drop in the function of a Reynolds number.
Consequently, the proposed numerical methodology resulted in a computational mesh that was orders of magnitude smaller when compared to the original 3D geometry. This allowed a wide range of geometrical modifications and flow conditions to be investigated. In the current study, the main focus is placed on reducing flow maldistribution, which is seen as a purely hydrodynamic phenomenon. Therefore, to minimize the influence of other factors, an incompressible and single-phase flow was considered.
It is important to note that in the previously motioned work [19], a two-phase flow with evaporation was modeled. A satisfactory comparison with an experiment confirmed that consistent geometrical reductions can be applied in the case of more complex flows.

Mathematical Model
The numerical model used in this work is based on an open source Computational Fluid Dynamics (CFD) numerical toolbox, OpenFOAM (Open Source Field Operation and Manipulation) [23,24]. This is a C++ object-oriented library for solving complex problems in CFD, and for carrying out structural analyses. The robustness and trustworthiness of this open source software has been proven multiple times in diverse and challenging applications. Among which, the works [25][26][27] deal with nonstandard and very complex flow and heat systems.
Numerical calculations were made using the simpleFoam solver available in OpenFOAM (release 1612+). This solver exploits finite volume discretization in conjunction with the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm for incompressible flows. According to OpenFOAM documentation [23], simpleFoam is a steady-state solver for incompressible flows with turbulence modeling.
The vector field of the flow velocity is found by solving the Navier-Stokes equation for a steady flow of incompressible fluid [28]: where u = (u, v, w) is a velocity vector, p is pressure and ν e f f is an effective kinematic viscosity calculated from selected transport and turbulence models.

Numerical Model and Transformation to 2D Geometry
PHEs usually consist of a large number of perforated plates. Assembled together, they form a complicated net of flow channels. The perforation of the plates usually exploit a chevron pattern, Figure 1. The very complex net of flow channels intensifies the turbulence and consequently enhances heat transfer. The current investigation is based on a real PHE device used in the experimental work of [12]. The geometrical details and measures are listed in Table 1. Figure 1. Comparison of real 3D PHE geometry and its 2D representation. For better visibility only a single pair of plates is shown. In the current studies a PHE assembled from 60 plates is considered. Substantial changes in the inlet and outlet headers geometry (collector pipes) and a lack of a chevron pattern in the 2D geometry can be observed. The geometrical complexity of the chevron plate results in a very complicated and dense numerical mesh. In the author's previous work [19], the minimum size of the 3D computational mesh for one flow channel (two plates), for a moderate Reynolds number, was shown to consist of ∼ 4.5 × 10 6 elements. A typical PHE has O(100) plates and its numerical model would require O(10 8 ) computational cells. Consequently, this would require an extremely lengthy amount of time for calculation and extensive hardware requirements.
Additionally, from an engineering perspective, the examination of a wide range of parameters is often necessary, with short computational time required. Consequently, full 3D computations can be prohibitively time and resource-consuming. This clearly motivates the need for a robust 2D reduction of the original 3D geometry.
An appropriate 2D reduction should preserve crucial dimensions, including: flow channel width W c , depth D (plate width) and the longitudinal height of the plates L HT . This assures the same heat transfer area and channel Reynolds number: where m in is the inletting mass flow, N c is the number of channels, D c is the hydraulic diameter of the channel and A c is a cross-sectional area and µ is viscosity. In the current study, a PHE with N c = 59 (60 plates) was considered. The PHE headers and flow channels were simplified to 2D but their cross-sectional areas were preserved to match the 3D geometry. Consequently, the inlet velocity consistency was satisfied. This requirement determined the height of the rectangular 2D inlet section to be: H in,out = 7.4 mm, see Figure 2a.
To reproduce the chevron flow pattern inside the PHE channels without a substantial increase of the mesh size, zero-thickness mixing baffles inside the 2D flow channels were introduced. They were placed alternately on both plates, as shown in Figure 2c. Applying the mixing baffles allowed a wavy flow pattern to be reproduced. This was previously shown to be a critical aspect of the 3D to 2D transition [19]. The amount of mixing baffles applied was determined by the analysis of the considered real 3D PHE geometry, and was estimated to be equal to 26. Another critical aspect of the 3D to 2D transformation, which has not been considered before, is related to the preservation of the pressure loss balance between the geometries. In the transformation from 3D to 2D geometry, there is a substantial change in the pressure loss balance. This is related to the different values for the hydraulic diameter of the header in the 3D geometry and its 2D representation. To even out the pressure loss balance, an originally designed approach has been proposed and is discussed in the next section. Figure 2. Sketch of the considered 2D model of PHE with additional geometrical modifications to assure similarity to the 3D geometry: (a) basic dimensions of 2D PHE; (b) placement of throttling baffles to preserve pressure balance between 2D and 3D geometries (blue); (c) mixing baffles to model the complexity of the 3D flow (red).

Pressure Balance Analysis and Throttling Baffles
An analysis of the proposed 2D geometry, when compared to its 3D original, shows a discrepancy in pressure loss balance across the PHE. Although the rectangular 2D header was modeled to preserve the inlet cross-section area, it did not conserve the same hydraulic diameter as with the circular 3D header. Additionally, due to a change in shape of the inlets to the PHE flow channels from circular (and localized in the middle of plates) to rectangular, there is a difference in the local pressure drop as well, see Figure 1. Consequently, there is a greater linear pressure drop in the header pipes (bottom and top collector pipe) but less local losses at the inlets and outlets of the PHE flow channels in the case of the 2D geometry, compared to the 3D geometry. A general idea to level these differences is to assure the same cumulative pressure losses in both PHE geometries: where ∆p 3D h and ∆p 2D h are the total pressure losses in the header of the 3D geometry and its 2D representation, respectively. ∆p 3D c and ∆p 2D c are the total pressure losses in the parallel channels of the PHE in the 3D geometry and its 2D representation, respectively. Further analysis assumes the same mass flow through each PHE channel (uniform flow distribution) and consequently, the same local and linear pressure loses in each individual channel.
It is justified to state that: The first inequality is obvious because the hydraulic diameter of the collector pipe of the 2D geometry is smaller than the hydraulic diameter of the corresponding 3D geometry, and can be proven by the Darcy-Weisbach formula. The second inequality is a hypothesis which will be checked numerically in the latter part of the paper.
Rearranging and expanding Equation (3) gives: c are the linear pressure losses in the header pipes of the 3D and 2D geometries, respectively, , ζ 3D in and ζ 2D in are the local pressure loss coefficients at the inlet to the PHE channels of the 3D geometry and its 2D representation, respectively. ζ 3D out and ζ 2D out are the local pressure loss coefficients at the outlet of the PHE channels of the 3D geometry and its 2D representation, respectively. L is the length of the collector pipe and λ is a friction factor To fulfill Equation (4) an additional pressure loss must be included to assure: The additional pressure losses: ζ 2D in,ex , ζ 2D out,ex and K 2D k,ex were included in the 2D geometry as properly designed local losses in the form of additional throttling baffles at the inlets and outlets of the PHE channels (Figure 2b).
The local loss coefficients are unknown and need to be evaluated numerically. For this purpose, two models of the PHE were created and each of them contained only two flow passages (four plates). The first geometry was based on the 3D PHE model and had a circular tube header (Figure 3a,b), and the second model was based on the 2D PHE and had a rectangular duct header (Figure 3c,d). To assure an accurate estimation of the local loss parameters, both models were created as three-dimensional, and each mesh contained over 30 million computational cells satisfying the mesh independence requirement.
The estimated local loss coefficients were based on the numerical results of the pressure and velocity fields measured on surfaces coloured in blue in Figure 3 and calculated using: ζ = 2∆p u 2 ρ . The obtained numerical values of the local loss coefficients for the considered range of Re c numbers are listed in Table 2. It can be seen, that the previously stated hypothesis: ∆p 3D c > ∆p 2D c is proven. Interestingly, the ratio of ζ 3D in, out \ ζ 3D in, out is independent from Re c and equals approximately 2.5 for the inlets and 2.0 for the outlets, despite the fact that the absolute value of the local loss coefficients have a tendency to decrease with Re c .  Based on Table 2 and Equation (5), it was possible to determine the values of the additional loss coefficients at the inlet and outlet ζ 2D in,ex and ζ 2D out,ex to fulfill Equation (4). This was incorporated in the numerical model with additional throttling baffles. Their length was calculated using values of in,ex and ζ 2D out,ex . To make this possible, it was necessary to determine the functional dependency of ζ coefficients from the throttling baffle size (the opening of an inlet and outlet of a PHE channel). In order to obtain these relations, a set of geometries with different baffle sizes was created. The influence of the bottom and top baffles were checked separately. The sizes of the inlet and outlet throttling baffle were changed gradually and the related pressure losses were calculated. The measurements of the pressure and velocity needed to calculate the local loss coefficients in the function of the channel opening were done on the same surfaces as before (Figure 3c,d). Their numerical values, expressed as percentages of the clearance of the inlet and outlet to the channel, are presented in Tables 3 and 4, respectively. The additional linear pressure loss in the header K 2D h,ex was calculated using the Darcy-Weisbach formula [29] and was directly included in the local losses coefficients. The final sizes of the throttling baffles placed at the inlet and outlet of the PHE channels are shown in Table 5. It can be observed, that the length of the throttling baffles gradually get smaller with the Re c . Table 5. Sizes of the additional throttling baffles in mm, which need to be placed at every inlet and outlet of the PHE channels to preserve the same pressure loss balance between the original 3D geometry and its 2D representation.

Variable
Re

Geometrical Modifications of PHE and Reduction of Flow Maldistribution
To decrease the intensity of the flow maldistribution and improve the performance of PHE, a number of different geometrical optimizations to the inlet header collector were examined. The main purpose of the proposed modifications was to force the fluid stream to enter the channels with a more uniformly distributed mass flow. The investigated geometrical changes are schematically shown in Figure 4 and summarized in Table 6. Their simplicity makes straightforward application possible in a real PHE. Consequently, all the configurations may be potentially applied in mass industrial production. Table 6. Considered geometrical optimization and explanation of the abbreviations used in this work, see Figure 4.

Geometry Description
org Basic 2D geometry of PHE without throttling and chevron baffles, Figure 2a.
org4 Developed 2D geometry of PHE with throttling baffles added at channel inlet and outlet to assure the same pressure balance in 2D and 3D geometries, Figure 2b. In the cases: 2str, 3str, 4str, 2str2, 3str2, 4str2 additional horizontal baffles (dividers) were introduced to the bottom header of the PHE to force a separation of the inletting fluid. Additionally, the influence of the placement of the starting point of the dividers was investigated, either from the first channel or from the beginning of the inlet header. These geometrical modifications were motivated by the expectation that the fluid would be redirected and delivered to different regions of the PHE. In the cases: TW-top, TW-alt vertical baffles were introduced instead of horizontal ones.
Their goal was to enhance turbulence and mixing in the bottom header. Theses baffles were placed in every fifth channel and covered half of the header cross-section. Figure 4. Schemes of the examined geometrical modifications. From left to right: 2str, 2str2, 3str, 3st2, 4str, 4str2, TW-top, TW-alt geometries, for a detailed explanation, see Table 6.

Results
In this work, flow maldistribution was evaluated based on a vertical component of the velocity vector in every channel of the PHE. The values of the velocity were read at a half height of the PHE. The mean velocity in every channel was calculated according to: where u y,i is the y component of velocity in the i-th channel.
To discuss and compare the flow distribution between different Re c numbers, the mass flow ratio coefficient was defined as: whereṁ i is the mass flow in the i-th channel andṁ ave is the theoretical average mass flow in the channel, assuming the perfect distribution of the flow. To evaluate and compare the results for different geometrical modifications, a degree of flow nonuniformity s was defined as an average of absolute deviation: where average mass flow ratio MFR ave is 1 by definition, and where s = 0 means there is an ideal flow distribution. Based on Equation (8), a degree of maldistribution parameter was defined as: where s w is the worst possible maldistribution calculated using Equation (8) with the assumption that the whole flow goes through the first channel of the PHE. The value of DM changes from 0 to 1, where 1 means a uniform mass flow in every channel (no maldistribution). Figure 5 shows the effect of Re c on the intensity of maldistribution. This compares the MFR coefficients for the considered 2D models: org, org4 and org5 for Re c = (50, 100, 250, 500, 1000). Note that the basic org geometry does not have mixing and throttling baffles and should be treated as reference results. The value of DM for these geometries are presented in Table 7.

Effect of Re c on Flow Maldistribution
The introduction of mixing and throttling baffles can be observed to have a crucial impact on the results. The differences between the org, org4 and org5 models are related to an inconsistent pressure drop balance as compared to the original 3D geometry (in the case of the org model). The pressure loss in the header of the org is too high in relation to the pressure drop in the channels, and the flow chooses "the easiest way" through the first few channels. Consequently, a higher flow maldistribution is observed. For every considered geometry, the intensity of the flow maldistribution grows with increasing Re c . This is particularly visible in the reference case org. It should be noted that these tendencies become substantially smaller in the more realistic model of the 2D PHE. For org5, this dependency is very weak, and can be seen in the bottom image of Figure 5. Nevertheless, the flow maldistribution is considerably strong in every case.     Tables 8-10 summarize the results of the considered geometrical modifications for the cases org, org4, org5, respectively. The geometries of org and org4 serve as a reference for the results of the full 2D PHE model for org5. Thanks to this, the influence of throttling and mixing baffles can be studied separately. Solutions containing additional horizontal baffles are observed to reduce the flow maldistribution significantly for the whole range of the considered Re c . It is not obvious what number of dividers is most preferable, but cases with a lower number of horizontal dividers tend to have more optimal outcomes. This is especially visible in the case of org5, which represents the full 2D model of the PHE. The lower DM value for a larger number of dividers can be related to the alternating levels of pressure balance in the PHE. Every additional divider creates a smaller hydraulic diameter in the lower collector pipe and, consequently, a larger pressure drop in the PHE header.
Opposite to this, in the case of vertical baffles (TW-top and TW-alt), the flow maldistribution becomes even greater when compared to the unmodified geometry of org5. The vertical obstacles create more turbulence in the header and the overall performance of the PHE seems to be less stable.
Consequently, this also causes a higher pressure drop in the PHE header and enhances the tendency of the flow to choose "the easiest way" through the first few channels.    Table 9. Values of the DM coefficient for the org4 geometry, Figure 2b.  In the case of vertical baffles, the flow is significantly blocked by the first baffle, which results in very similar PHE performance for both considered solutions with vertical baffles, compared with the last row in Figure 8.

Conclusions
The presented work is focused on developing an efficient and fast numerical methodology for comprehensive research related to PHEs. The robust transformation of the real 3D geometry of a PHE into a 2D representation was a key feature in order to significantly reduce the time for computation, as well as to provide trustworthy results. The proposed geometrical transformation reduced the number of computational mesh cells by two orders of magnitude. A key feature of the 2D transformation was the introduction of mixing and throttling baffles. The number and placement of the mixing baffles were directly related to the original 3D geometry and their task was to mimic a wavy flow pattern caused by the chevron pattern. The necessity of the throttling baffles was justified to assure the same pressure balance between the 3D and 2D geometries. Their size was calculated based on the assumption that a cumulative pressure drop should be preserved between the original 3D geometry and its 2D representation. The developed 2D model made it possible to perform a significant number of calculations, and to explore a wide range of parameter space, Re c = (50, 100, 250, 500, 1000), and a number of various geometrical modifications introduced to the PHE. The main idea of the proposed geometrical changes was to assure its simplicity and ease of production. Consequently, a number of cases introducing additional horizontal and vertical dividers to the header of the PHE (bottom collector pipe) was considered. An analysis of the results provided a potentially meaningful insight into more preferable geometrical modifications leading to a more uniform flow distribution in the PHE channels.
In the case of a typical PHE, a reasonable solution is to introduce horizontal dividers to its bottom collector. This divides the inlet flow into several streams and directly delivers it into farther sections of the PHE. More optimal solutions are with a smaller number of horizontal dividers. Every additional divider creates a smaller hydraulic diameter and, consequently, a larger pressure drop in the PHE header. In the case of introducing vertical baffles to the PHE header, flow maldistribution becomes even greater than with a PHE without any baffles. The vertical obstacles intensify the turbulence and the overall performance of the PHE seems to be less stable. Consequently, this causes a higher pressure drop in the PHE header too. The above conclusions should also hold in the case of two-phase flows or flows with a phase change.
The proposed methodology of the 2D transformation can be potentially seen as universal. It should be possible to adopt the developed methodology to different types of PHE. Different perforation patterns of the PHE plates could be potentially modeled with different distribution, sizes and angles of mixing baffles. In general, they should mimic the original plate pattern "as closely as possible". In the case of a lower Reynolds numbers, this should provide trustworthy results because the flow in a real PHE is approximately two-dimensional. In the author's previous work [19] the flow in-between the plates was shown to remain two-dimensional up to a channel Reynolds number of 250 for V-chevron plates. Consequently, the main limitation of the proposed 2D methodology can be disclosed in the case of a higher Reynolds numbers, because the flow in a real PHE would become fully three-dimensional. Assuming that a flow with a higher Reynolds number would mostly influence the overall pressure balance in the PHE, the 2D model can be adjusted to the 3D original by assuring the same pressure balance. This can be achieved using a similar methodology to the one presented in this work. It would require full 3D calculations, but only for one flow channel (single pair of plates).
It is worth noting that the proposed methodology is purely geometrical and independent from a mathematical model. It is ready to be used for a multi-phase flow or a flow with a phase change. In such cases, the size of the throttling baffles should be determined using an adequate flow model, but the general methodology would be similar to the one presented in the current studies.