NUMERICAL INVESTIGATION OF TUBESIDE MALDISTRIBUTION IN SHELL-AND-TUBE HEAT EXCHANGERS

Shell-and-tube heat exchangers (STHE) are widely used in the process and energy industry. Maldistribution and the often resulting fouling in these STHE cause additional energy consumption and lower production throughput. Increased average wall shear stress, compared to that of the maldistributed case inside the tubes is expected to mitigate fouling. This can be achieved by a uniform distribution into the tubes. Field and laboratory data suggest that common crude oil fouling is profoundly mitigated above 10 Pa and is significantly reduced above a wall shear stress of 15 Pa [1]. Many STHE already have lower design shear stresses than those mentioned. Therefore, if maldistribution takes place, tubes with less flow velocity will have even more fouling. To investigate tubeside flow maldistribution, a parametric STHE model is studied with computational fluid dynamics (CFD). At first, a comparison between the standard k-𝜖 -model and the new standard SST-model is performed to check if SST could provide improved simulation results. Afterward, a range of geometrical parameters will be investigated to find influencing quantities of maldistribution. The resulting velocity distributions are visualized and evaluated by using different statistical approaches. At least, a sensitivity analysis will be done to show how each parameter influences the tubeside flow distribution in STHE.


INTRODUCTION
Shell-and-tube heat exchangers are used in all areas of energy and process engineering processes to transfer heat flows from one fluid to another. According to [2] STHE have a market share of more than 60 % of all heat exchangers used. Therefore, every optimization of the heat transfer performance has the potential to reach a very large group of users worldwide.
Operating experiences typically report fouling on the tube and shell side of a large number of these devices (up to 90 % of all heat exchangers [3]), which have to be cleaned with high energy and financial effort [4,5]. Furthermore, the systems work less efficiently, the operating costs increase and there are additional emissions of greenhouse gases.
Fouling on the tubeside is caused by fluid properties, unfavorable design and disadvantageous operating parameters, and by maldistributed flow through the tube bundle, which in turn is favored by unfavorable design.
Uniform flow through the tube bundle represents the thermo-fluid-dynamic ideal condition and is preferred to minimize the fouling tendency of several fluids and fluid systems. Standardized design guidelines for STHE, for example, according to the standards of the Tubular Exchanger Manufacturers Association (TEMA) [6] or the VDI Heat Atlas [3], generally assume an even flow, which occurs merely randomly and is not the norm. This leads to a higher fouling tendency and general performance losses. Design standards used today are essentially based on empirical studies, empirical values and assumptions, whereby some of them are decades old. They are thus suitable for designing a range of apparatus as wide as possible but cannot address specifics of individual materials or process properties.
Reviewing these standards with modern tools like CFD (computational fluid dynamics) brings new insights and fundamental knowledge, which can help increase the efficiency of these assets and make them less prone to fouling.
In 1988, Mueller and Chiou [7] published a paper summarizing all their work on maldistribution up to that time. They found that an axially mounted inlet nozzle leads to a jet impingement into the tube sheet direction. Therefore, it was assumed that a radially mounted inlet nozzle could avoid the jet impingement and lead to a uniform flow distribution. That was investigated by several studies, such as Kim et al. [8], Mohammadi et al. [9], and Dorau et al. [10]. All concluded that a radial mounted inlet nozzle could not prevent maldistribution. Mueller and Chiou further concluded that most heat exchanger applications would not notice the performance reduction according to maldistribution as it results in minor losses (< 5 %) when the number of transfer units (NTU) is below 10. However, if NTU exceeds the NTU = 10 criteria, losses up to 15 % can occur [7].
Bremhorst and Brennan [11] performed CFD studies to investigate the influence of flow phenomena in 4-pass-STHE and their impact on tube inlet wear, so the metal loss at the tube entrances of the tube bundle. One part of their study compared the standard k--model against a Large-Eddy-Simulation (LES) to predict the flow through the first pass of the STHE. It was found that the LES provides more accurate results. On the other hand, the authors concluded that the improved accuracy achieved was insufficient to justify the computational effort as the flow distribution of both models showed very similar tendencies. They continued with the k--model.
By using the SST-model (Shear Stress Transport) instead of the k--model, an improvement in simulation accuracy is expected as it allows a much better resolution of the flow boundary layer near the wall. This can lead to a more accurate prediction of the flow behavior inside a STHE and especially in the tube bundle.
In 1977 Gotoda and Izumi [12] investigated the influence of the attachment angle of the inlet nozzle on the flow distribution in the tube bundle of a ship condenser. In their experiments, they found that the maximum deviation from the mean velocity could be estimated by using equation (1).
Here describes the flow velocity with max for maximum deviated velocity and 0 for design flow velocity. While is the dimensionless total loss coefficient and the different cross-section area (A) subscripts stand for 1 = inlet pipe, 2 = tube sheet and 3 = tube bundle.
Using the equation, they predicted that an area ratio AR ~ 1 (equation (2)) should produce a uniform flow in the tube bundle. The same was found as one main influence factor for maldistribution throughout a CFD parameter study with several refinery STHEs [13]. Said et al. [14] investigated how orifices and nozzles (tube nozzles) at the inlet of each tube of a cross-flow heat exchanger (HX) affect the velocity distribution. They found that orifices can reduce maldistribution by up to 12 times, as it works like a second tube sheet or a perforated plate mentioned here [7], while the overall pressure drop was increased by about 8 %. Tube nozzles, on the other hand, could reduce maldistribution only about 7.5 times but with an additional decrease of pressure drop of almost 10 %. As this study aimed to decrease both maldistribution and pressure drop, the tube nozzle approach is tested for STHE applications.

Fouling through Maldistribution
Longstaff and Palen [1] found in field and laboratory experiments threshold criteria where asphaltenic fouling can be mitigated at a wall shear stress of approximately 10 Pa and is strongly minimized when it exceeds 15 Pa.
In addition, Joshi et al. [15] suggested the utilization of a fouling model without the removal. This model was developed based on the investigation of several crude oils in a double-pipe heat exchanger test plant. By neglecting the proposed aging and surface roughness terms, the following simplified power law approach between wall shear stress and fouling rate can be drawn: where A and B are oil specific constants. This approach allows the description of the underlying fouling phenomena when maldistribution occurs. Similar shear stress dependent fouling behaviors can be found for crystallization and biofouling [3,16].
Assuming an average design shear stress of 10 Pa into the tubes, the positive and negative deviation of flow velocities will lead to several problems. By summarizing the results of [1] and [15] it can be concluded, that with the appearance of maldistribution, tubes with shear stress below the 10 Pa threshold (see Figure 2). This will lead to an exponential increase of fouling respectively, the lower the flow velocity into the tubes is. On the other hand, tubes with high flow velocities will have only slightly reduced fouling and might be harmed otherwise, for example, by erosion and corrosion phenomena [11,17].

Performance loss through Maldistribution
STHE are typically designed with the assumption of uniform flow distribution in the tube bundle. Several analytical studies [7,10,18,19] pointed out, that this assumption might be good enough for a broad range of heat exchanger applications but can also bring performance losses up to 15% right from the start of operation. tubes are x % slower than the average velocity [10] Dorau et al. [10] compared the performance loss under clean and fouled conditions after one year of operation. They used a simplified heat exchanger model with artificial flow distributions in the range from uniform to up to 50 % of tubes with 50 % higher velocities (vmax = 3.0 m/s) than the theoretical average tube velocity v0 = 2.0 m/s) and vice versa (vmin = 1.0 m/s). In this study the flow distribution had only an impact on the tube outlet temperature and thus on the logarithmic mean temperature difference. Figure 3 shows this performance loss through maldistribution. While ̇ decreases by about 8 %, ̇ was decreased by 12 % over the same range of maldistribution. Further comparing the difference of heat duty loss of ̇ against ̇, it can be concluded that the fouled case had lost more heat duty over the range of maldistribution.
In summary, previous investigations show that flow maldistribution in the tube bundle has a negative influence on heat transfer and fouling behavior, especially in tubes with low flow velocities. At the same time, there is a risk of mechanical damage for pipes with very high flow velocities, for example, by increased erosion or cavitation due to sudden pressure drop at tube entrances [11], which in total means considerable costs for plant operators. Therefore, reasons for maldistribution in STHE have been investigated by CFD, which will be described in the following section.

Numerical modeling
Basic modeling parameters of the CFD setup and information about the evaluation workflows (heat maps, bar charts, and validation of the shortened setup) can be found in [13]. As before, only the tubeside inlet nozzle, inlet header, and the first meter of the tube bundle have been modeled to decrease the computational effort. The validity of this approach was already verified in the previous work [13]. The current study was conducted with fluid parameters of an industrial oil as introduced in Said et al. [14], since the results of both studies will be compared later. The simulation was set up isothermally, which is typical for such numerical flow studies. The density was 645.6 kg/m3, and the dynamic viscosity was 2.32·10-3 Pas. Both values are in the range of crude oil at bulk temperatures (50 °C to 250 °C) in refinery pre-heat-trains (PHT) [20,21].
In addition, the current work was conducted using a STHE-model constructed according to typical design standards [3] with a number of tubes of n = 54. As an extension to the previous CFD model [13] this CAD (computer aided design) geometry (see Figure 1 and Figure 4) was parametrized to check all specific variables uncoupled from each other. About 100 different parameter sets (simulation runs) have been investigated. As before, ANSYS® CFX was used to perform the CFD simulations. The mesh was proofed by a mesh independence study at the basic geometry first and for both RANS models. The element count was about 50 million volume elements for each set (meshing setup similar to [13]). However, at first, the mentioned comparison of different numerical models is presented. Comparison of − vs. SST Bremhorst and Brennan [11] investigated if and how the simulation results will be improved by using LES modeling instead of the standard k-model for flow simulation inside of the tubes of a heat exchanger.
Instead of using LES, it was investigated here if the usage of the standard k--model would still be sufficient for the aims of the study compared to the upcoming new standard SST model. With a sufficiently small meshing of the areas close to the wall (so-called prism layer), the SST model enables a direct calculation of the velocities present in the boundary layer instead of only approximating them with a logarithmic wall function, as is the case with k-. The criterion for "sufficiently small" in this case is a dimensionless wall distance + < 5 (~1 at best).
On the other hand, using SST means much more computational effort for solving the same problem and leads to additional difficulties in getting a wellconverged solution. With the higher resolution of the SST model, it was impossible to get a simple steady-state solution, as turbulences were resolved on a lower scale. In consequence, the targeted residuals are not achievable or only with difficulty. When this happens, there might be two approaches. The first approach would be to reduce the mesh density, which would destroy the improvement of the model's more exact physics. The authors decided to use a secondary approach. Instead of a steadystate simulation, a transient approach has been set up. The resulting values (mainly velocity and pressure) were averaged over 60 s after the simulation to get comparable values with the steadystate k--solution. For the transient simulation a time step of 0.15 s (Δ ) was used. Figure 5 shows the comparison of the different flow velocity curves over time of the transient solution (blue), the averaged solution (red), and the steady-state k-solution (yellow) of one randomly chosen tube of the tube bundle. It can be seen that the averaged transient solution deviates only about 2 to 3 % from the steady-state results (in each tube), which was considered negligible. The same applies to the pressure drop values. Therefore, it was decided, as already described by Bremhorst and Brennan, that the standard k--model is still sufficient for such a broad parameter study.

Parametric investigation of header design
In this section, the most promising results of the parameter study will be shown. Single results are plotted against the uniformity index 2 (4), as introduced in [13] respectively [22], and the resulting pressure drop calculated from CFD.
At the end a summary of all results will then be compared, as a sensitivity analysis, with the linear correlation coefficient (5) as introduced in [13].
The velocity distribution of the basic design ( Figure 6 and Figure 7) already shows a good distribution ( 2 = 0.9806). This could later be described by the relatively large LK and aS values, which positively influence the distribution (see Figure 8). Typical values for 2 vary between 0.9 (bad) and 0.99 (almost perfect). Four test series showed that a short inlet header with low aS is unfavorable for the uniform flow through the tube bundle. Increasing the length of the inlet chamber has a positive effect, although in most cases, the impact is even higher if the nozzle is also moved to the back. In some cases, small local maxima for 2 occurred as a function of aS when LK/aS = 2. The global maxima, at constant LK, were achieved in every case where the nozzle was placed as far as possible from the tube sheet. Figure 8 shows the comparison of the velocity magnitude on the symmetry plane of two parameter sets (a) (LK = 500 mm, aS = 400 mm) and (b) (LK = 300 mm, aS = 150 mm) in detail. In the area between the inlet nozzle and the tube sheet, the velocity varies highly over the height of the inlet header. As the distance from the nozzle plane increases, the velocity profile becomes more uniform. This phenomenon occurs over the complete cross-section of the header. Since aS is larger in (a) than in (b), the flow velocity immediately upstream of the tube sheet is much more uniform. The values of 2 combined with the graphically displayed simulation results show that the velocity profile in front of the tube bundle influences the flow through the individual tubes. The more uniform it is, the more similar the single flow velocities are within the tubes. The flow into the header from the normal direction to the tube bundle also turned out to be an unfavorable inlet nozzle configuration (Figure 9), as the distribution will be focused only on a certain amount of tubes directly opposite the inlet nozzle. This configuration leads to 2 of 0.95, which is one of the worst values in this study. Izumi and Gotoda [12] mentioned that an reduction of AR (eq. 2) would lead to a more uniform flow distribution across the tube bundle, this was also found in [13]. Having an area ratio AR≤1 could be achieved by decreasing A3 (cross-sectional area of the tube bundle), by increasing the cross-sectional area of the inlet nozzle A1 or by both methods in parallel.
A reduction of A3 has a positive effect in terms of 2 and thus reduces the disadvantage of an increased pressure drop somewhat ( Figure 10). Especially, since Δ would increase even with other common methods for fouling reduction. The problem with changing the cross-sectional area of the tube bundle is that this also affects the thermal behavior of the heat exchanger. This complicates the design of the apparatus, for which many criteria already have to be taken into account. An increase of the nozzle diameter ( Figure 11) on the other hand would be a much simpler approach. It leads to lower inlet velocities, which in turn results in lower Reynolds number and therefore to lower turbulence and pressure drop into nozzle and header. Furthermore, the lower inlet velocity reduces the jet impingement on the pass partition plate. Which results in smaller flow velocity deviations near the plate and the flow has more time to align into tube bundle direction. In conclusion this means that increasing the nozzle diameter would be the preferable method to reach the AR≤1 criteria.
Due to the knowledge that the pressure drop is high at the tube entrance a combination of the discussed methods (A1 ↑ or A3 ↓) may further enhance the advantages of both individual measures.

Figure 11: Uniformity index and pressure drop over varied inlet nozzle diameter (DS)
Said et al. [14] used nozzles at the tube inlet of each tube (see Figure 12) of a smaller heat exchanger with nine tubes and achieved very good results in terms of flow uniformity and pressure drop. This approach was investigated here as well. The inner diameter of the tube bundle tubes remains at Di = 20 mm, while the diameter at the tube sheet is conical increased up to DTN = 24 mm. The length of the tube nozzle LTN is the same as in the investigations by Said et al. [14] (LTN = 20 mm) and is equal to the inner diameter of the tubes.
In order to assess whether adding tube nozzles to the heat exchanger tubes represents a reasonable change in geometry, the results of the "Tube nozzle" test ( Figure 13) were compared with the test series with varied inner tube diameters, in particular with the Di = 20 mm test setup ( Figure 10). The inner tube diameter is the same in both cases. Thus, the change in Δ and 2 by the installed tube nozzles could be investigated in detail. The improvement in Δ seems to be larger than the deterioration in 2 . Thus, this simulation results suggest that by using tube nozzles, it is possible to mitigate the negative influence of smaller inner tube diameters Di on the pressure drop without reducing their positive effect on the uniformity index to the same extent. This single investigation carried out, does not allow to determine if this is the case in general. Since the attachment of tube nozzles at the beginning of the tube bundle would be costly regarding manufacturing technology, the flow properties would have to be significantly improved to justify this effort. In the investigation carried out, a positive influence of the tube nozzles could be observed, but not to the same extend as reported in [14].
Investigations undertaken by Dorau et al. [13] showed that the velocity of the inflow into the header only has a minor influence on the flow distribution in the tube bundle. Therefore, the impact of a load change (regarding the flow velocity) mostly tends to affect the spread of the tube-discrete maximum or minimum velocities. In turn, this means that an optimization of the in-and outflow conditions of the inlet header simultaneously makes the heat exchanger less sensitive to load changes. The flow uniformity is thus maintained at low to medium load deviations from the design conditions. This was confirmed by the current study (see Figure  14). It can be seen that increasing the mass flow rate will decrease the uniformity index only slightly while increasing the pressure drop significantly.  Figure 15 shows the sensitivity analysis where the impact of several parameters on the uniformity index is compared by the correlation coefficient rxy (5).

Figure 15: Bar chart with correlation coefficient over certain varied parameters
The influence of the chamber length LK and the distance of the nozzle to the tube sheet aS is difficult to distinguish. In conclusion, larger LK and aS lead to higher the uniformity indices.
Varying (the rotation angle of the inlet nozzle) and bS (the lateral point of nozzle attachment at the header) have provided better values of 2 only in isolated cases. In most cases, they worsened the uniform flow through the tube bundle.
In contrast, the variation of inlet nozzle and tube diameters (presented as area here) and the associated values like v and Re show mostly negative or no significant impacts. The tube bundle cross-section (A3) shows the greatest statistical influence on the distribution, but when it comes to economic considerations there might be a tradeoff with the highly increased pressure drop.
As mentioned before, an adjustment of the area ratio AR (eq. 2) would be the best way to decrease the maldistribution. It could be shown that an increased AR would result in an increased maldistribution.
In conclusion, this work confirms the results of the previous study [13]. It can be said that this is a further step in investigating maldistribution on the tubeside of STHE. A heat exchanger model (CFD/ analytical) with added heat transfer and/or dynamic fouling build up would provide deeper insights into the underlying mechanisms.

SUMMARY
Using preliminary work and key literature references, the influence of tubeside maldistribution on STHE performance and fouling behavior was reviewed. Both phenomena reduce the originally intended (design) performance and increase maintenance effort and costs. With this in mind, numerical simulations were carried out to identify the basic geometrical factors influencing maldistribution.
A comparison of the two RANS models k-and SST should serve as a basis to determine whether the increased numerical effort of SST compared to kleads to a significant improvement of the simulation results. This could not be proven for the present application. So the standard k--model was used for the following study.
Subsequently, a wide-ranging parameter study was carried out. It turned out that many of the investigated parameters only have a small influence on the tubeside maldistribution, while a few parameters could bring a larger improvement against flow maldistribution.
In particular, the distance of the inlet nozzle from tube sheet aS, the length of the inlet chamber LK and the area ratio AR should be mentioned, which, when applied correctly, produce a significant increase in the uniform distribution with a simultaneous reduction in pressure drop.

NOMENCLATURE
A material specific constant, m 2 K/W/d A1 cross-section area of inlet nozzle, m 2 A2 tube sheet area, m 2 A3 cross-section area of tube bundle, m 2 a distance, m B material specific constant, dimensionless b width/offset, m D diameter, m L length, m ̇ mass flow rate, kg/s n number of tubes, dimensionless R thermal resistance, m 2 K/W rxy correlation coefficient, dimensionless t time, s v flow velocity, m/s x statistical variable, dimensionless y statistical variable, dimensionless y + dimensionless wall distance, dimensionless uniformity index, dimensionless total loss coefficient, dimensionless dynamic viscosity, Pas density, kg/m 3 rotation angle, ° shear stress, Pa