Nonlinear Finite Element Analysis-Based Flow Distribution and Heat Transfer Model †

: A new strategy for fast, approximate analyses of ﬂuid ﬂow and heat transfer is presented. It is based on Finite Element Analysis (FEA) and is intended for large yet structurally fairly simple heat transfer equipment commonly used in process and power industries (e.g., cross-ﬂow tube bundle heat exchangers), which can be described using sets of interconnected 1-D meshes. The underlying steady-state model couples an FEA-based (linear) predictor step with a nonlinear corrector step, which results in the ability to handle both laminar and turbulent ﬂows. There are no limitations in terms of the allowed temperature range other than those potentially stemming from the usage of ﬂuid physical property computer libraries. Since the ﬂuid ﬂow submodel has already been discussed in the referenced conference paper, the present article focuses on the prediction of the tube side and the shell side temperature ﬁelds. A simple cross-ﬂow tube bundle heat exchanger from the literature and a heat recovery hot water boiler in an existing combined heat and power plant, for which stream data are available from its operator, are evaluated to assess the performance of the model. To gain further insight, the results obtained using the model for the heat recovery hot water boiler are also compared to the values yielded by an industry-standard heat transfer equipment design software package. Although the presented strategy is still a “work in progress” and requires thorough validation, the results obtained thus far suggest it may be a promising research direction.


Introduction
During design, operation, and troubleshooting of various process and power equipment-containing tube bundles, it often is important to know the velocity and temperature fields in both the tube and the shell sides. These are obtained predominantly using Computational Fluid Dynamics (CFD) models and, therefore, articles covering a wide range of such applications are available. For example, Wei et al. [1] discussed a coupled CFD-Lagrange multipliers optimization method for flow distribution adjustments to prevent freezing of power generation natural draft cooling systems during winter operation. Chien et al. [2], on the other hand, presented a coupled CFD-surrogate-based optimization of flow distribution in a heat exchanger inlet header. Zhou et al. [3] focused on CFD investigation and optimization of a compact heat exchanger comprising a single row of tubes, and Łopata et al. [4] published an article covering the experimental investigation of flow distribution in a similar cross-flow heat exchanger, but with a tube bank consisting of elliptical tubes. CFD evaluation and optimization of solar collectors, commonly also using a single row of risers, were discussed, for instance, by García-Guendulain et The original model discussed in [26] assumed adiabatic flow, that is, no heat transfer was allowed on the walls of the parallel flow channels in the distribution system. The model was shown in the same article to provide data with relative errors of at most 4% compared to detailed transient CFD simulations even in the case of highly turbulent flows. Such accuracy was achievable due to the relative simplicity of the flow systems for which the respective model has been intended (e.g., tube bundles in heat recovery steam generators). The overall conclusion, therefore, was that, in terms of application in preliminary analyses of selected heat transfer equipment or for shape optimization of the mentioned equipment, the model was suitable for engineering practice.
Because of the nature of the model, its performance in case of laminar flow was a priori expected to be acceptable. Although several tests were carried out earlier even with low total mass flow rates to make sure this really was the case, no example was given in [26]. To remedy this, let us mention, for instance, one of the test flow distribution systems (see Figure 1) used in the original article and the respective laminar flow distribution data and relative errors. For convenience's sake, parameters of the flow system are listed in Table 1. The obtained mass flow rates are compared in Figure 2a, while Figure 2b shows the corresponding relative errors. It can be seen that the error values generally were in a ±1% band with only two of them being at around 1.2%. Relative errors obtained using other test flow systems were of similar magnitudes. Thus, one could conclude that, in the case of laminar flow, the accuracy was even better than when the flow was highly turbulent.  Table 1. Figure 1. Schematic of the flow distribution system from Table 1. Table 1. Parameters of the flow distribution system used to obtain the laminar flow-related data shown in Figure 1.

Parameter Value
headers (W × H × L) 55 × 55 × 280 mm tube bundle 5 rows with 10 tubes each, tube layout: 60 • tubes straight, inner diameter: 10 mm, length: 2000 mm fluid water, 0.5 kg s −1 , 300 K flow regime laminar, average tube Re ≈ 1500 Figure 1. Schematic of the flow distribution system from Table 1.  Table 1 using the model based on Finite Element Analysis (FEA) discussed in [26] and a transient Computational Fluid Dynamics (CFD) simulation. Average tube Reynolds number was ca. 1500. (b) The corresponding relative tube mass flow rate errors (FEA vs. CFD simulation). Tube numbers correspond to Figure 1. For the details regarding the CFD model, the reader is kindly referred to [27].

Inclusion of Heat Transfer into the Model
The main shortcoming of the original, flow-only version of the Finite Element Analysis (FEA)based model was its inability to properly evaluate tube bundles in which heat transfer could not be neglected. Given the intended purpose of the model (that is, usage in engineering practice), this functionality had to be implemented.
Please note that heat transfer was not, strictly speaking, evaluated using FEA. However, the temperature fields in the tube side and the shell side were still obtained using a system of linear equations generated as shown in the following text, and this system was then solved in the usual manner. It was assumed that the temperature profile between two end nodes of an edge was continuous and was given by the mean temperatures on control volume cross-sections, which were perpendicular to the corresponding edge. The iterative solver then worked similarly to any other segregated solver. First, the fluid flow (FEA-based) submodel was solved under the assumption of a constant temperature field. Next, the heat transfer submodel was solved under the assumption of a constant velocity field. This was followed by the update of the fluid physical properties and other necessary post-iteration tasks, and then the fluid flow submodel was solved again. This iterative procedure was repeated until convergence was reached.  Table 1 using the model based on Finite Element Analysis (FEA) discussed in [26] and a transient Computational Fluid Dynamics (CFD) simulation. Average tube Reynolds number was ca. 1500. (b) The corresponding relative tube mass flow rate errors (FEA vs. CFD simulation). Tube numbers correspond to Figure 1. For the details regarding the CFD model, the reader is kindly referred to [27].

Inclusion of Heat Transfer into the Model
The main shortcoming of the original, flow-only version of the Finite Element Analysis (FEA)-based model was its inability to properly evaluate tube bundles in which heat transfer could not be neglected. Given the intended purpose of the model (that is, usage in engineering practice), this functionality had to be implemented.
Please note that heat transfer was not, strictly speaking, evaluated using FEA. However, the temperature fields in the tube side and the shell side were still obtained using a system of linear equations generated as shown in the following text, and this system was then solved in the usual manner. It was assumed that the temperature profile between two end nodes of an edge was continuous and was given by the mean temperatures on control volume cross-sections, which were perpendicular to the corresponding edge. The iterative solver then worked similarly to any other segregated solver. First, the fluid flow (FEA-based) submodel was solved under the assumption of a constant temperature field. Next, the heat transfer submodel was solved under the assumption of a constant velocity field. This was followed by the update of the fluid physical properties and other necessary post-iteration tasks, and then the fluid flow submodel was solved again. This iterative procedure was repeated until convergence was reached.
Even though the heat transfer submodel was not using FEA, the corresponding mesh on which the temperature field was calculated can be constructed in a similar manner. In the fluid flow mesh, the field to be calculated was described by total pressures in the nodes. The temperature field can be described analogously with the difference being that every edge had its own temperature in the node. Figure 3 shows the two meshes and the differences between meshes. Even though the heat transfer submodel was not using FEA, the corresponding mesh on which the temperature field was calculated can be constructed in a similar manner. In the fluid flow mesh, the field to be calculated was described by total pressures in the nodes. The temperature field can be described analogously with the difference being that every edge had its own temperature in the node. Figure 3 shows the two meshes and the differences between meshes. As mentioned before, the temperatures were calculated using a set of linear equations. From Figure 3 it is clear that a flow system consisting of n edges will feature 2n node temperatures and, therefore, 2n linear equations were required. There were three classes of temperature-related equations that were used in the model:  Flow mixing and splitting,  Heat transfer through channel walls, and  Boundary conditions.

Flow Mixing and Splitting
When, in an arbitrary mesh node, streams q1, q2, …, qm are mixed into a single stream j, we can write Here, ̇ denotes the mass flow rate of the qth stream, cp,q the specific heat capacity, and Tj and Tq the corresponding stream temperatures. Each specific heat capacity should be taken as the mean value obtained for the corresponding temperature range [Tj, Tq]. If, on the other hand, a single stream j is split into streams r1, r2, …, rn, the outflow temperature is the same for all these streams, and the respective n equations are In some systems, there may be blind edges with zero mass flow rate. The temperatures in the nodes of these edges are calculated as if the edges were of the regular type featuring outflow (see also the schematic in Figure 4).  As mentioned before, the temperatures were calculated using a set of linear equations. From Figure 3 it is clear that a flow system consisting of n edges will feature 2n node temperatures and, therefore, 2n linear equations were required. There were three classes of temperature-related equations that were used in the model:

•
Flow mixing and splitting, • Heat transfer through channel walls, and • Boundary conditions.

Flow Mixing and Splitting
When, in an arbitrary mesh node, streams q 1 , q 2 , . . . , q m are mixed into a single stream j, we can write Here, . m q denotes the mass flow rate of the qth stream, c p,q the specific heat capacity, and T j and T q the corresponding stream temperatures. Each specific heat capacity should be taken as the mean value obtained for the corresponding temperature range [T j , T q ].
If, on the other hand, a single stream j is split into streams r 1 , r 2 , . . . , r n , the outflow temperature is the same for all these streams, and the respective n equations are T r = T j , r ∈ {r 1 , r 2 , . . . , r n }. ( In some systems, there may be blind edges with zero mass flow rate. The temperatures in the nodes of these edges are calculated as if the edges were of the regular type featuring outflow (see also the schematic in Figure 4). Even though the heat transfer submodel was not using FEA, the corresponding mesh on which the temperature field was calculated can be constructed in a similar manner. In the fluid flow mesh, the field to be calculated was described by total pressures in the nodes. The temperature field can be described analogously with the difference being that every edge had its own temperature in the node. Figure 3 shows the two meshes and the differences between meshes. As mentioned before, the temperatures were calculated using a set of linear equations. From Figure 3 it is clear that a flow system consisting of n edges will feature 2n node temperatures and, therefore, 2n linear equations were required. There were three classes of temperature-related equations that were used in the model:  Flow mixing and splitting,  Heat transfer through channel walls, and  Boundary conditions.

Flow Mixing and Splitting
When, in an arbitrary mesh node, streams q1, q2, …, qm are mixed into a single stream j, we can write ∑̇p , ( − ) Here, ̇ denotes the mass flow rate of the qth stream, cp,q the specific heat capacity, and Tj and Tq the corresponding stream temperatures. Each specific heat capacity should be taken as the mean value obtained for the corresponding temperature range [Tj, Tq]. If, on the other hand, a single stream j is split into streams r1, r2, …, rn, the outflow temperature is the same for all these streams, and the respective n equations are = , ∈ { 1 , 2 , … , }.
In some systems, there may be blind edges with zero mass flow rate. The temperatures in the nodes of these edges are calculated as if the edges were of the regular type featuring outflow (see also the schematic in Figure 4).  In a general case with streams q 1 , q 2 , . . . , q m being mixed and then split into streams r 1 , r 2 , . . . , r n , one will get one Equation (1) governing the resulting outflow temperature T j and (n − 1) Equation (2), Energies 2020, 13, 1664 6 of 20 that is, (n − 1) identities for the remaining outflow temperatures. The total number of equations governing the mixing/splitting in the node will, therefore, be equal to number of outflow streams.

Heat Transfer through Channel Walls
Let us have two meshes representing the tube and the shell sides of a heat exchanger and focus on an arbitrary pair of adjacent control volumes representing a portion of the tube side (i.e., a tube segment) and the enclosing portion of the shell side (see Figure 5). Let . m t , c p,t , T t,1 , and T t,2 denote the tube side mass flow rate, mean specific heat capacity at constant pressure, and inlet and outlet temperatures and . m s , c p,s , T s,1 , and T s,2 the corresponding shell-side quantities.
Energies 2020, 12, x FOR PEER REVIEW 6 of 20 In a general case with streams q1, q2, …, qm being mixed and then split into streams r1, r2, …, rn, one will get one Equation (1) governing the resulting outflow temperature Tj and (n − 1) Equation (2), that is, (n − 1) identities for the remaining outflow temperatures. The total number of equations governing the mixing/splitting in the node will, therefore, be equal to number of outflow streams.

Heat Transfer through Channel Walls
Let us have two meshes representing the tube and the shell sides of a heat exchanger and focus on an arbitrary pair of adjacent control volumes representing a portion of the tube side (i.e., a tube segment) and the enclosing portion of the shell side (see Figure 5). Let ̇t, cp,t, Tt,1, and Tt,2 denote the tube side mass flow rate, mean specific heat capacity at constant pressure, and inlet and outlet temperatures and ̇s, cp,s, Ts,1, and Ts,2 the corresponding shell-side quantities. Should, e.g., the hot fluid be placed in the shell side, then the overall heat balance could be written as ̇t p,t ( t,2 − t,1 ) =̇s p,s ( s,1 − s,2 ). ( Let us for a moment assume that the temperature of the fluid in the shell-side control volume is constant and is equal to the shell-side inlet temperature, Ts,1. Let us also assume that the tube-side inlet temperature, Tt,1, is known. Additionally, let L denote the length of the tube-side mesh edge and U the cumulative overall heat transfer coefficient. The heat flux for a small portion of this edge having the length dx can then be expressed as This can be modified, rearranged, and written in integral form, which yields From Equation (6), we immediately see that the temperature at the end of the edge can be obtained using Equation (7) must be linearized for it to be used in a matrix solver. This is done in a straightforward manner by taking Should, e.g., the hot fluid be placed in the shell side, then the overall heat balance could be written as Let us for a moment assume that the temperature of the fluid in the shell-side control volume is constant and is equal to the shell-side inlet temperature, T s,1 . Let us also assume that the tube-side inlet temperature, T t,1 , is known. Additionally, let L denote the length of the tube-side mesh edge and U the cumulative overall heat transfer coefficient. The heat flux for a small portion of this edge having the length dx can then be expressed as This can be modified, rearranged, and written in integral form, which yields From Equation (6), we immediately see that the temperature at the end of the edge can be obtained using Energies 2020, 13, 1664 7 of 20 Equation (7) must be linearized for it to be used in a matrix solver. This is done in a straightforward manner by taking where the necessary cumulative overall heat transfer coefficient, U, is computed from the tube-side and shell-side heat transfer coefficients. These, in turn, are calculated using equations from literature (e.g., [28] in the case of plain tubes or [29] if the tubes are finned) depending on the actual bundle geometry and flow conditions. Additional information regarding validation of the commonly used empirical equations for estimation of heat transfer coefficient in the case of plain and serrated fins can be found, for instance, in [30]. One could also use the correlations from [31], which have been obtained via the machine learning technique. If U-shaped or helical fins of various types are employed, then the two-part article by Hofmann and Heimo [32,33] can be recommended to the reader. The final, linearized equation for a single edge, therefore, is Considering the shell-side outlet temperature, for cross-flow with T s,1 = const. on the entire inlet face of the control volume, we have Just as before, this can be modified and rearranged to yield The mean shell-side outlet temperature then is which corresponds to the shell-side outlet temperature obtained using the respective set of linear equations.
One could also simplify the model even further by using a one-dimensional shell-side mesh (i.e., a mesh such that each cross-section of the shell along the general flow direction is spanned by just one cell). With a row of n tubes being present in a specific shell-side cell, Equation (3) while Equation (7), still necessary for each of the n tubes, would remain almost identical: There also is a special case of no heat transfer, which can be treated similarly. The necessary equation can easily be obtained by setting the heat transfer coefficient in Equation (14) to zero, which results in the equation being reduced to the equality between the temperatures in the end nodes of an edge. This is important because the number of linear equations describing heat transfer is always constant no matter if heat transfer occurs or not.

Boundary Conditions and the Complete System of Linear Equations
Up to this point, every equation was simply describing some relationship between the node temperatures. For a steady state problem to be completely specified, some temperatures must be known. However, let us first analyze the number of equations available thus far. For a fluid flow system with n edges and m inflow streams, there are 2n unknown temperatures. We can get n equations from the heat transfer. The following n − m equations can be obtained from stream mixing in the nodes. The remaining m equations must be provided via boundary conditions, i.e., inlet temperatures must be specified for each of the inflow stream (other arrangements may be possible in selected cases). When there are multiple fluid flow systems connected by heat transfer equations, the number of available equations remains the same.

Coupling of the Flow Distribution and Heat Transfer Submodels
Each major iteration of the FEA-based solver consists of two steps. The first step is a fixed temperature field analogy of the adiabatic model (as described in [26]; robustness of the model can be improved by carrying out this first step repeatedly until the respective residual falls below a predefined threshold). New estimates of the temperature fields for both the tube and the shell sides are then obtained in the second step. Here, the necessary values of C U are updated using edge mass flow rates from the first step and the corresponding new estimates of cumulative overall heat transfer coefficients. To solve the respective combined linear system for the tube-side and the shell-side temperatures, one boundary temperature must be provided for each stream (for instance, at the inlet of each tube in the bundle and for each inlet cell in the discretized shell side). The resulting temperature matrix could look, for example, like the one in Figure 6. Even though linear systems represented by such matrices can be solved quite easily, it is obvious from the figure that implementation of a reordering algorithm would be necessary should one want to improve performance via preconditioning in case of much larger linear systems. temperatures must be specified for each of the inflow stream (other arrangements may be possible in selected cases). When there are multiple fluid flow systems connected by heat transfer equations, the number of available equations remains the same.

Coupling of the Flow Distribution and Heat Transfer Submodels
Each major iteration of the FEA-based solver consists of two steps. The first step is a fixed temperature field analogy of the adiabatic model (as described in [26]; robustness of the model can be improved by carrying out this first step repeatedly until the respective residual falls below a predefined threshold). New estimates of the temperature fields for both the tube and the shell sides are then obtained in the second step. Here, the necessary values of CU are updated using edge mass flow rates from the first step and the corresponding new estimates of cumulative overall heat transfer coefficients. To solve the respective combined linear system for the tube-side and the shell-side temperatures, one boundary temperature must be provided for each stream (for instance, at the inlet of each tube in the bundle and for each inlet cell in the discretized shell side). The resulting temperature matrix could look, for example, like the one in Figure 6. Even though linear systems represented by such matrices can be solved quite easily, it is obvious from the figure that implementation of a reordering algorithm would be necessary should one want to improve performance via preconditioning in case of much larger linear systems. Figure 6. Non-zero values in the sparse temperature field matrix used in the second step of the FEA solver. Please note that, for the sake of clarity, only a small matrix with the rank slightly below 800 is shown here, which originates from a flow distribution system with two bundles consisting of just four tubes each.
As the convergence criterion, the fluid flow submodel used the scaled norm of the difference between the solutions from the predictor step and the corrector step. The corresponding scaled residual limit was 10 −5 . In the case of the heat transfer submodel, if we denote the original linear system Ax = b, then the scaled residual norm is computed from b − Ax just before the heat transfer submodel is solved. (If it were done after the respective solution process, the norm would be equal to zero.) The same residual limit, that is, 10 −5 , was used here.
All physical property data (mean specific heat capacity, dynamic viscosity, etc.) are always taken Figure 6. Non-zero values in the sparse temperature field matrix used in the second step of the FEA solver. Please note that, for the sake of clarity, only a small matrix with the rank slightly below 800 is shown here, which originates from a flow distribution system with two bundles consisting of just four tubes each.
As the convergence criterion, the fluid flow submodel used the scaled norm of the difference between the solutions from the predictor step and the corrector step. The corresponding scaled residual Energies 2020, 13, 1664 9 of 20 limit was 10 −5 . In the case of the heat transfer submodel, if we denote the original linear system Ax = b, then the scaled residual norm is computed from b − Ax just before the heat transfer submodel is solved. (If it were done after the respective solution process, the norm would be equal to zero.) The same residual limit, that is, 10 −5 , was used here.
All physical property data (mean specific heat capacity, dynamic viscosity, etc.) are always taken for the current conditions from the IAPWS [34] or the CoolProp [35] physical property libraries, or, in special cases (e.g., flue gas), are computed using various interpolation functions or tabulated data depending on the actual compositions. Thermal properties of the tube and fin materials are always obtained using tabulated data from literature (for example, from the technical standard [36]).

Shell-Side Pressure Drop
Similarly to heat transfer coefficients, pressure drop in the shell side cross-flow tube bundle is calculated via well-known empirical equations from, e.g., [37]. The actual formula to be used depends on the bundle geometry, possible presence of fins, etc.

The Developed Computer Code
The computer code was developed in Python [38] and utilized NumPy [39] to carry out the necessary matrix computations. The Visualization Toolkit (VTK) [40] and meshio [41] libraries were used to export solution data to Kitware ParaView [42] for visualization. Although no graphical user interface (GUI) is available yet, the authors plan to add it in the future, for example, via the Django web framework [43]. Please note that the code is not publicly available.
Apart from the inclusion of heat transfer, many additional improvements of the code have been made since the publication of the initial article discussing the FEA-based model [26]. The most important one probably is parallelization of the mass flow rate corrector step (please see [26] for details). As the correction algorithm was carried out independently for each mesh edge, a set of asynchronous workers was created using the standard Python multiprocessing library, and the correction tasks were processed in batches on all available CPU cores. This then results in much shorter computational times. Parallelization of the internal matrix computations, however, was not implemented because, given the numbers of elements in the simplified meshes, the matrices were rather small. In other words, the CPU time saved by parallel solution would be wasted on auxiliary operations needed to split the task to multiple cores, thus rendering the net improvement either negligible or even negative.

Results
In this section, two test cases are discussed to demonstrate the capabilities of the present version of the developed model. First, a simple cross-flow tube bundle heat exchanger from the literature is evaluated in Section 3.1. Section 3.2 then focuses on a heat recovery hot water boiler in an existing combined heat and power plant, for which stream data have been provided by its operator.

Simple Cross-Flow Tube Bundle Heat Exchanger from the Literature
The example discussed here was taken from [44] and involves an air-to-water heat exchanger from Figure 7. Its parameters are listed in Table 2 together with the data obtained using the present model and HTRI Xchanger Suite 8.0.1 [45]. The computational time needed by the present model to automatically create the necessary meshes, reach in 46 major iterations the results mentioned below, and export all the solution data into Kitware ParaView for visualization purposes was ca. 15 s on an average desktop computer with the Intel Core i-5 2500K CPU. The ranks of the matrices used in the model were ca. 800 and ca. 1600 in case of fluid flow and heat transfer, respectively.  In order to minimize the number of sources of discrepancies, the necessary heat transfer coefficients were calculated by the present model using the equations mentioned in the literature [44]. The results should, therefore, have been identical, yet they were not. The reason for the difference became clear once one noticed that, in [44], the iterative computation was stopped while the difference between the hot and the cold heat duties was still relatively large. In fact, should one carry out the heat balance for the data from the literature, one would get the actual water heat duty of 150.3 kW (as listed among the results) while for air the duty would be 149.9 kW. There may also be another reason for the discrepancies in the data, namely the fact that, in [44], the computation was done using average temperatures, average fluid physical properties, etc. for the entire tube side and shell side.
Considering the differences between the values provided by the present model and the data yielded by HTRI Xchanger Suite, these most likely stemmed from the software package using much more accurate equations for obtaining the heat transfer coefficients. It, therefore, is entirely possible that with different equations the results obtained using the present model would be much closer to the data from HTRI Xchanger Suite. In any case, this supports the notion that any such tool or model can only be as good as the equations internally utilized by it.
To demonstrate the level of detail of the solution data provided by the developed software, the water and air temperatures were exported (together with the respective geometry, which used only plain tubes to improve clarity) into Kitware ParaView. The resulting combined plot is shown in Figure 8.  Table 2. Parameters of the air-to-water heat exchanger (for the remaining construction data, etc., please see [44]) and the corresponding results obtained using the present model and HTRI Xchanger Suite ("HTRI XS").

Parameter
Literature [ In order to minimize the number of sources of discrepancies, the necessary heat transfer coefficients were calculated by the present model using the equations mentioned in the literature [44]. The results should, therefore, have been identical, yet they were not. The reason for the difference became clear once one noticed that, in [44], the iterative computation was stopped while the difference between the hot and the cold heat duties was still relatively large. In fact, should one carry out the heat balance for the data from the literature, one would get the actual water heat duty of 150.3 kW (as listed among the results) while for air the duty would be 149.9 kW. There may also be another reason for the discrepancies in the data, namely the fact that, in [44], the computation was done using average temperatures, average fluid physical properties, etc. for the entire tube side and shell side.
Considering the differences between the values provided by the present model and the data yielded by HTRI Xchanger Suite, these most likely stemmed from the software package using much more accurate equations for obtaining the heat transfer coefficients. It, therefore, is entirely possible that with different equations the results obtained using the present model would be much closer to the data from HTRI Xchanger Suite. In any case, this supports the notion that any such tool or model can only be as good as the equations internally utilized by it.
To demonstrate the level of detail of the solution data provided by the developed software, the water and air temperatures were exported (together with the respective geometry, which used only Energies 2020, 13, 1664 11 of 20 plain tubes to improve clarity) into Kitware ParaView. The resulting combined plot is shown in Figure 8.

Heat Recovery Hot Water Boiler in an Existing Plant
A heat recovery hot water boiler with nominal thermal power of 53.3 MW in an existing combined heat and power (CHP) plant was selected as the second test case. The boiler contained two counter-flow tube bundles which were mounted in the vertical portion of the flue gas duct (i.e., the tubes are horizontal, see Figure 9). Both bundles consisted of several passes, and each pass was composed of four staggered tube rows with 48 tubes per row. In the bottom bundle, the first pass was unfinned, the second pass used plain round fins, and serrated fins were utilized in the third pass. The top bundle contained solely tubes enhanced with serrated fins. The built-up area of the heated portion of each bundle was ca. 7.6 × 4.0 m. All stream-related data presented in this article were obtained by the operator of the boiler in the course of a guarantee test.

Heat Recovery Hot Water Boiler in an Existing Plant
A heat recovery hot water boiler with nominal thermal power of 53.3 MW in an existing combined heat and power (CHP) plant was selected as the second test case. The boiler contained two counter-flow tube bundles which were mounted in the vertical portion of the flue gas duct (i.e., the tubes are horizontal, see Figure 9). Both bundles consisted of several passes, and each pass was composed of four staggered tube rows with 48 tubes per row. In the bottom bundle, the first pass was unfinned, the second pass used plain round fins, and serrated fins were utilized in the third pass. The top bundle contained solely tubes enhanced with serrated fins. The built-up area of the heated portion of each bundle was ca. 7.6 × 4.0 m. All stream-related data presented in this article were obtained by the operator of the boiler in the course of a guarantee test.
The boiler was driven by flue gas exiting from a gas turbine. Because the temperature field (see the measurement array shown in Figure 9) was almost uniform, the corresponding boundary condition was specified in both the models discussed further as a constant. The outlet temperature of flue gas was estimated by the operator because the respective quantity had not been measured during the guarantee test. All the necessary data are summarized in Table 3. tubes are horizontal, see Figure 9). Both bundles consisted of several passes, and each pass was composed of four staggered tube rows with 48 tubes per row. In the bottom bundle, the first pass was unfinned, the second pass used plain round fins, and serrated fins were utilized in the third pass. The top bundle contained solely tubes enhanced with serrated fins. The built-up area of the heated portion of each bundle was ca. 7.6 × 4.0 m. All stream-related data presented in this article were obtained by the operator of the boiler in the course of a guarantee test.  The simulation carried out using the developed computer code mentioned in Section 2.3 included both water distribution in the tube side and heat transfer between the flue gas and water. To assess the accuracy of the predicted temperatures, the boiler was also analyzed using an industry-standard tool, namely, HTRI Xchanger Suite. Please note that, with respect to the requests of the manufacturer of the boiler and the operator of the CHP plant, no other data regarding the apparatus can be explicitly specified in this article. For the same reason, neither the HTRI Xchanger Suite case files nor the simplified 3-D CFD model of the flue gas duct discussed in the following sections can be made available.

Simulation in HTRI Xchanger Suite
Compared to a full-scale CFD simulation of the boiler, which would rarely be done in the case of equipment of such size, the actual computational time required by HTRI Xchanger Suite was negligible (units of seconds). Unlike CFD, however, the software generally focuses on the thermal side of the apparatus design, i.e., its primary goal is to provide as accurate stream temperatures as possible while the flows in both the tube and shell sides are assumed to be uniformly distributed (unless the user specifies the distribution explicitly). Moreover, one cannot directly set tube inner and outer surface roughnesses, which may significantly influence the predicted pressure drops.
Results obtained using the discussed software package are listed in Table 4 together with the data provided by the operator of the heat recovery hot water boiler. The table also mentions the absolute and relative errors. From these, one can see that the predictions of both the outlet temperatures and the tube-side pressure drop were quite accurate. The predicted shell-side pressure drop, however, was markedly lower than the measured value. It also was a notable disadvantage that no detailed information was given by HTRI Xchanger Suite regarding the actual flow distributions in the bundles and the shell. Table 4. Results obtained using HTRI Xchanger Suite and the corresponding errors compared to the data from the operator of the boiler. It is apparent that the temperatures and the tube-side pressure drop were predicted quite accurately, but the shell-side pressure drop was markedly below the measured value.

Assessment of the Shell-Side Flow Behavior
In order to verify whether the assumption of uniform velocity distribution over the flue gas duct cross-section was appropriate in the FEA-based computation discussed in the next section, a simplified 3-D CFD model of the duct was created in ANSYS Fluent [46]. Parameters of the model are listed in Table 5. To keep computational demand at a reasonable level, the bundles were replaced by porous zones. Additionally, the entire duct was split into several parts, and hexahedral cells were used whenever possible to further lower the cell count while maintaining acceptable mesh quality. After the necessary mesh adaptation, the final cell count was ca. 3.3 M (see the y + histogram in Figure 10). As for cell equiangle skewness, only 1231 cells (ca. 0.04% of the total number of cells) featured skewness greater than 0.6, of which 1227 fell between 0.6 and 0.7. The obtained contour plot of velocity magnitude just below the bottom bundle is then shown in Figure 11. This indicates that, although the velocity distribution was not entirely uniform (see also the pathlines in Figure 12), the non-uniformity was still at a reasonable level, which should not lead to significant inaccuracy in the calculated overall heat transfer rate.  In order to verify whether the assumption of uniform velocity distribution over the flue gas duct cross-section was appropriate in the FEA-based computation discussed in the next section, a simplified 3-D CFD model of the duct was created in ANSYS Fluent [46]. Parameters of the model are listed in Table 5. To keep computational demand at a reasonable level, the bundles were replaced by porous zones. Additionally, the entire duct was split into several parts, and hexahedral cells were used whenever possible to further lower the cell count while maintaining acceptable mesh quality. After the necessary mesh adaptation, the final cell count was ca. 3.3 M (see the y + histogram in Figure  10). As for cell equiangle skewness, only 1231 cells (ca. 0.04% of the total number of cells) featured skewness greater than 0.6, of which 1227 fell between 0.6 and 0.7. The obtained contour plot of velocity magnitude just below the bottom bundle is then shown in Figure 11. This indicates that, although the velocity distribution was not entirely uniform (see also the pathlines in Figure 12), the non-uniformity was still at a reasonable level, which should not lead to significant inaccuracy in the calculated overall heat transfer rate.   Table 5. There were 275 cells with y + ≤ 30 (all of these but 12 were between 25 and 30, ca. 0.04% of the total number of wall-adjacent cells) and no cells with y + > 300.  Table 5. There were 275 cells with y + ≤ 30 (all of these but 12 were between 25 and 30, ca. 0.04% of the total number of wall-adjacent cells) and no cells with y + > 300.
Energies 2020, 12, x FOR PEER REVIEW 14 of 20 Figure 11. Contour plot of velocity magnitude just below the bottom bundle, which was obtained using the CFD model from Table 5. The velocity distribution was not uniform (due to the flue gas duct being bent) and fluctuated in time.

Figure 12.
Pathlines in the flue gas duct, which were colored by velocity magnitude. One can see that the flow behind the bend is pushed to one side as indicated in Figure 11 and that there are relatively large recirculation zones present.

Present Model
The data yielded by the present, FEA-based model are listed in Table 6 together with the values provided by the operator of the boiler. Here, the accuracy was slightly lower than that of HTRI Xchanger Suite, but it still was acceptable. The computational time needed to automatically create the meshes, reached in 39 major iterations the solution, and export the necessary solution data into Kitware ParaView for visualization purposes was ca. 240 s on the same average desktop computer used in the previous example. The ranks of the matrices were ca. 9000 and ca. 18,500 in case of fluid flow and heat transfer, respectively. Table 6. Results obtained using the present model and the corresponding errors compared to the data from the operator of the boiler. The accuracy was not as good as in the case of HTRI Xchanger Suite, but, in terms of fast, approximate analyses of process and power equipment, it was sufficient.  Figure 11. Contour plot of velocity magnitude just below the bottom bundle, which was obtained using the CFD model from Table 5. The velocity distribution was not uniform (due to the flue gas duct being bent) and fluctuated in time.

Parameter
Energies 2020, 12, x FOR PEER REVIEW 14 of 20 Figure 11. Contour plot of velocity magnitude just below the bottom bundle, which was obtained using the CFD model from Table 5. The velocity distribution was not uniform (due to the flue gas duct being bent) and fluctuated in time.

Figure 12.
Pathlines in the flue gas duct, which were colored by velocity magnitude. One can see that the flow behind the bend is pushed to one side as indicated in Figure 11 and that there are relatively large recirculation zones present.

Present Model
The data yielded by the present, FEA-based model are listed in Table 6 together with the values provided by the operator of the boiler. Here, the accuracy was slightly lower than that of HTRI Xchanger Suite, but it still was acceptable. The computational time needed to automatically create the meshes, reached in 39 major iterations the solution, and export the necessary solution data into Kitware ParaView for visualization purposes was ca. 240 s on the same average desktop computer used in the previous example. The ranks of the matrices were ca. 9000 and ca. 18,500 in case of fluid flow and heat transfer, respectively. Table 6. Results obtained using the present model and the corresponding errors compared to the data from the operator of the boiler. The accuracy was not as good as in the case of HTRI Xchanger Suite, but, in terms of fast, approximate analyses of process and power equipment, it was sufficient.  Figure 12. Pathlines in the flue gas duct, which were colored by velocity magnitude. One can see that the flow behind the bend is pushed to one side as indicated in Figure 11 and that there are relatively large recirculation zones present.

Present Model
The data yielded by the present, FEA-based model are listed in Table 6 together with the values provided by the operator of the boiler. Here, the accuracy was slightly lower than that of HTRI Xchanger Suite, but it still was acceptable. The computational time needed to automatically create the meshes, reached in 39 major iterations the solution, and export the necessary solution data into Kitware ParaView for visualization purposes was ca. 240 s on the same average desktop computer used in the previous example. The ranks of the matrices were ca. 9000 and ca. 18,500 in case of fluid flow and heat transfer, respectively.
As mentioned before, the present model uses a one-dimensional mesh to represent the shell side. This means that the predicted temperature distribution was, too, only one-dimensional, while flow distribution across the shell-side cannot be predicted at all. In the tube side, on the other hand, the mass flow rate was known for each individual tube, and the predicted temperature distribution was spatially as fine or as coarse as the utilized tube bundle mesh. Figure 13 shows the predicted shell-side temperature distribution along the portion of the flue gas duct enclosing the two bundles, which was obtained using the present model, and the corresponding temperatures provided by the operator of the boiler and yielded by HTRI Xchanger Suite. The temperature curve yielded by the model matches the  Table 6. Results obtained using the present model and the corresponding errors compared to the data from the operator of the boiler. The accuracy was not as good as in the case of HTRI Xchanger Suite, but, in terms of fast, approximate analyses of process and power equipment, it was sufficient.

Parameter
Value As mentioned before, the present model uses a one-dimensional mesh to represent the shell side. This means that the predicted temperature distribution was, too, only one-dimensional, while flow distribution across the shell-side cannot be predicted at all. In the tube side, on the other hand, the mass flow rate was known for each individual tube, and the predicted temperature distribution was spatially as fine or as coarse as the utilized tube bundle mesh. Figure 13 shows the predicted shellside temperature distribution along the portion of the flue gas duct enclosing the two bundles, which was obtained using the present model, and the corresponding temperatures provided by the operator of the boiler and yielded by HTRI Xchanger Suite. The temperature curve yielded by the model matches the point values reasonably well with the discrepancies being most likely caused by the usage of different equations for the calculation of the necessary heat transfer coefficients. Figure 13. Comparison of the temperature distribution in the portion of the flue gas duct enclosing the two bundles, which was obtained using the present model, and the corresponding temperatures provided by the operator of the boiler and yielded by HTRI Xchanger Suite ("HTRI XS"). Vertical distance along the channel corresponds to the distance from the point denoted "y = 0 m" in Figure 9.
Similarly as before, a combined plot of water and flue gas temperatures was generated using Kitware ParaView. This is shown in Figure 14. The obtained tube mass flow rates were distributed quite uniformly in both bundles. The actual relative standard deviations from uniform flow distribution computed using Figure 13. Comparison of the temperature distribution in the portion of the flue gas duct enclosing the two bundles, which was obtained using the present model, and the corresponding temperatures provided by the operator of the boiler and yielded by HTRI Xchanger Suite ("HTRI XS"). Vertical distance along the channel corresponds to the distance from the point denoted "y = 0 m" in Figure 9.
Similarly as before, a combined plot of water and flue gas temperatures was generated using Kitware ParaView. This is shown in Figure 14.
The obtained tube mass flow rates were distributed quite uniformly in both bundles. The actual relative standard deviations from uniform flow distribution computed using where . m id denotes the ideal mass flow rate through one tube of the bundle, n the number of tubes therein, and . m i the mass flow rate through the ith tube, were 0.3% and 0.4% in case of the top bundle and the bottom bundle, respectively. Figure 13. Comparison of the temperature distribution in the portion of the flue gas duct enclosing the two bundles, which was obtained using the present model, and the corresponding temperatures provided by the operator of the boiler and yielded by HTRI Xchanger Suite ("HTRI XS"). Vertical distance along the channel corresponds to the distance from the point denoted "y = 0 m" in Figure 9.
Similarly as before, a combined plot of water and flue gas temperatures was generated using Kitware ParaView. This is shown in Figure 14.

Discussion
The results yielded by the present model and by HTRI Xchanger Suite for the air-to-water heat exchanger from the literature have shown that even though an apparatus can be decomposed into parts for which correlations or calculation procedures for heat transfer coefficients and pressure losses may exist, successfully applying them may not be straightforward. The main reason is that such procedures require local fluid and material properties, and these generally depend on the temperature and pressure, that is, quantities which the designer is trying to calculate. This is where the present model steps in. The data have also highlighted the facts that the accuracy of any model depends to a large degree on the quality of equations utilized for the calculation of the various coefficients and that further research in this area is, therefore, necessary. Additionally, one can draw the conclusion that using averaged values for the entire tube side and shell side can lead to notable differences. In the case of the outlet temperatures in this air-to-water heat exchanger, it was up to ca. ±7%.
As for the heat recovery hot water boiler, the most important stream parameters are listed in Table 7. It can be seen that the tube-side outlet temperature and pressure drop were better predicted by HTRI Xchanger Suite, while the accuracy of the shell-side outlet temperature and pressure drop predictions was higher in the case of the present model. Overall, the accuracy of the model was deemed acceptable concerning its intended purpose. Table 7. Summary of the main results obtained using the two discussed approaches alongside the data provided by the operator of the boiler. of the tubes in the bundle and the subsequent mechanical failures). Another advantage is that in the case of cross-flow tube bundles (the primary target application of the present model), HTRI Xchanger Suite assumes that mixing of tube streams occurs after each pass, even if this may not actually be true. Such a simplification may increase convergence, but it also may diminish local effects and, therefore, introduce errors into the data.

Conclusions
The point of this research was not to develop a replacement to the universally applicable, commercial heat transfer equipment design packages, such as HTRI Xchanger Suite. On the contrary, the goal was to create a simplified model for heat exchangers representable using sets of interconnected 1-D meshes, which are typically used in high-temperature (i.e., heat recovery) industrial applications and are prone to suffer from operating problems. The model, once finished, should be fast, yet accurate enough, and should provide estimates of not only the flow distribution in the bundle and the tube-and shell-side temperature fields but also the resulting mechanical stress field in the bundle caused by uneven thermal loading. In other words, the aim was to have a supplementary tool which would enable the designer to evaluate in advance, and without any significant effort or time spent, the thermal-hydraulic behavior of the mentioned heat exchangers as well as the likelihood of them suffering mechanical failures under the design operating conditions. In this regard, the FEA-based modelling approach seems to be promising, but a lot of work, as well as a thorough validation, are still needed before it is ready for production use.
The comparison with HTRI was mentioned in this paper solely to present the current capabilities of the model in terms of heat transfer prediction. The discrepancy (or at least a part of it) was very likely caused by the heat transfer coefficients and the hydraulic resistance coefficients being calculated differently. However, the model was designed in such a way that it can easily incorporate any standard method for calculating these coefficients for the simple 1-D mesh elements. More complex flow systems can then be built from these simple elements and the solution strategy remains the same, which ensures scalability of the model.
Considering the fact that shell-side flow velocity fields commonly are not uniform, one of the possible future improvements of the FEA-based model could lie in strictly using a grid of cells in the shell side (in the plane perpendicular to the general direction of flow) instead of only one cell. Another enhancement, which the authors plan to implement, is to make it possible to interface the current computer code with other simulation codes. This would enable the evaluation of flow systems in which some parts are more complex and, therefore, not directly compatible with the simplified mesh elements used by the FEA-based model.  Reynolds number, -T temperature, • C U cumulative overall heat transfer coefficient, W K −1 x length coordinate, m y + dimensionless wall distance, -Greek Symbols: δ relative standard deviation form uniform flow distribution, % ε rate of dissipation of turbulent kinetic energy, m 2 s −3 Subscripts: 1 at the inlet of the control volume 2 at the outlet of the control volume i, j, q, r summation indices id ideal value s in the shell side t in the tube side