Document Type : Research Article
Authors
a K. N. Toosi University of Technology, Faculty of Mechanical Engineering, DOS Computational Research Lab, Tehran, Iran
Abstract
Keywords
Main Subjects
Declarations
Funding
This work was partially supported by the National Iranian Gas Company, NIGC [grant number 950406].
Conflict of Interest Statement
The authors declare that they have no affiliation, association, or financial involvement with any organization or commercial entity with a direct or indirect financial interest or conflict with the subject matter or research presented in the manuscript.
This work was partially supported by the National Iranian Gas Company, NIGC [grant number 950406] but included no contractual or implied restriction on utilization or publication of the data and/or review of the data before publication.
Availability of data and material
Numerical data resulting in this paper is available via (ShafieiAlamooti & Ashrafizadeh, 2020).
Article Highlights
Acknowledgments
Funding: This work was supported by the National Iranian Gas Company, NIGC [grant number 950406].
Nomenclature 

Area (m^{2}) 
A 

Base coefficient in the reconstruction () 


Effective width of HE (m) 
B 

Correlation matrix () 
C, c_{pq} 

Distance in the Re vector space () 
d_{s} 

Relative information content () 
E_{1n} 

Drift in the KRG method () 
F 

Fin height (m) 
H 

Thermal conduction coefficient (W/mK) 
K 

Effective length (m) 
L 

Number of snapshots () 
M 

Number of nodes, Degree of freedom () 
N 

Pressure (N/m^{2}) 
P 

Heat transfer (W) 
Q 

Reynolds number () 
Re 

Actual range in the KRG method () 
r_{a} 

Fin pitch (m) 
S 

Number of the nearest points in the IDW method () 
S 

Temperature (K) 
T 

Matrix of the reconstructed field () 


Snapshot vector () 


Snapshot vector component () 


Components of velocity (m/s) 
u,v, w 

Specific internal energy (j/kg) 
u 

Velocity vector (m/s) 


Volume defined in a general vector space () 
v 

Matrix of weighting coefficients for the snapshot () 
W 

Weighting coefficient of a base for the reconstruction of a snapshot () 
w 

Fin inside width, Geometrical coordinate (m) 
x 

Fin inside height, Geometrical coordinate (m) 
y 

Departure in the KRG method () 
Z 

Geometrical coordinate (m) 
z 

Greek letters 

Semivariance () 


Deviation from CFD result () 
Δ 

Fin and separating walls thickness (m) 
δ 

Relative deviation percentage () 


Dimensionless temperature () 


Actual range percentage () 


Eigenvalue () 


Exponent in the IDW method () 


Density (kg/m^{3}) 
ρ 

Variance () 


Base vector () 


Matrix containing all base vectors () 


Subscripts 

Longitudinal 
L, long. 

Heat exchange 
Q 

Temperature 
T 

Superscripts 

Aluminium 
Al 

Average value 
ave 

Cold stream 
cold 

Desired reconstruction 
D 

Hot stream 
hot 

Intermediate stream 
int. 

Of a known coefficient 
s 

Dimensionless 
* 

Mean value 


Fluctuation value in the affine subspace 


MultiStream Heat Exchangers (MSHEs) have widespread applications, especially in cryogenics (Das & Ghosh, 2012) and gas liquefaction (Ikealumba & Wu, 2014). They are the key elements of cold boxes in LNG refineries (Zargoushi et al., 2020). Their high compactness and adaptability for energy integration also lead to a further complicated heat transfer pattern in comparison with the conventional twostream ones; thus, there is a high demand for costeffective analysis methods for them (Das & Ghosh, 2012). Mostly, substantial simplifications are inevitable in the presented methods. MSHEs are initially analyzed by Haseler (HASELER, 1983) and Paffenbarger (Paffenbarger, 1990). They used some correlations for the convection coefficients and estimated the temperatures using the fin efficiency concept. Chiba (Chiba, 2010) presented a method for conjugate thermal analysis of fluids and solids, i.e., the Graetz problem, in an MSHE with laminar flow. Krishna et al. (Krishna et al., 2013) showed the effect and importance of longitudinal heat conduction in them. Phasechange heat transfer processes of an MSHE are modeled by Niroomand et al. (Niroomand et al., 2020). They used the available correlations for convection coefficients and predicted the temperature distribution in the heat exchanger body. The proposed method is also utilized to investigate the effect of property variation on the flow maldistribution (Niroomand et al., 2021). CFD approach has not been so conventional for multistream heat exchangers because of its expensive computational requirements (Huang et al., 2014). Zargoushi et al. (Zargoushi et al., 2020) performed such a numerical study and considered only one flow condition for modeling an LNG refinery cold box. It has two thermal communications, and boiling and condensation are taken into account. Haider et al. (Haider et al., 2020) simplified a largescale MSHE by using established correlations for the interaction between fluid flow and fin structure and presented a transient numerical simulation via OpenFOAM v6. They combined this simulation with a genetic algorithm optimization routine to obtain a feasible design with high efficiency, low pressure drop, and low unit weight (Haider et al., 2021). Shafieialamooti and Ashrafizadeh (Shafieialamooti & Ashrafizadeh, 2021) did a phenomenological analysis on a platefin MSHE with three thermal communications. They investigated the interaction between property variation and thermal parasitic effects for several flow scenarios. The considered parasitic effects are collectively entitled Internal Heat Leak mechanisms.
As shown, the high computational requirements of the detailed numerical analysis of MSHEs make it necessary to accept substantial simplifications and use costeffective techniques. Soft computing methods like Neural Networks and Genetic algorithms are frequently used in conventional twostream heat exchanger modeling and optimizations. Aasi and Mishra (Aasi & Mishra, 2021) used an artificial neural network modeling to predict the thermal and hydraulic parameters such as Colburn and friction factor for a platefin heat exchanger. Yin and Ooka (Yin & Ooka, 2015) simulated a passage of a platefin heat exchanger to find correlations for the convection coefficient along the fins. They used it in the optimization procedure via a genetic algorithm. Proper Orthogonal Decomposition (POD) is an order reduction method effectively used for the study of thermofluid systems, such as data center server cooling systems (Ghosh & Joshi, 2014), automotive HVAC systems (Christ & Sattelmayer, 2018), thermoacoustic heat engine (Selimefendigil & Öztop, 2014), and PhaseChange Material (PCM) (Selimefendigil & Öztop, 2021). POD analysis has potential applications for MSHEs, which have not been considered yet. There are very few POD approaches for heat exchangers in the literature. Du et al. (Du et al., 2013) studied an aircooled condenser using POD. They investigated the effects of wind velocity on the performance of the condenser. Wang et al. (Wang et al., 2017) examined the variation of three design parameters for a finned tube heat exchanger. Both of these two studies only modeled the flow outside the tubes.
POD utilizes some available solutions called snapshot vectors to approximate the intended field at an acceptable cost. A PODaided simulation tries to extract a set of orthogonal bases, 's first. Then the field is modeled as a linear combination of them, , hereafter simply called the reconstruction. The main steps of the procedure are (Holmes et al., 1998; Sirovich, 1987):
For a thermofluid problem, POD analysis is strongly affected by the selection of the configuration of orthogonal bases: the scalarvalued approach (Brenner et al., 2009; Rowley et al., 2001) or the vectorvalued approach (Brenner et al., 2009; Rowley et al., 2001). The former is a straightforward choice but requires more weighting coefficients. The latter has a more complicated procedure for the computation of the bases but requires fewer weighting coefficients.
The intended unknown field is indicated by parameters like a particular time, geometrical dimensions, values at boundaries, etc. The snapshots must be collected for various values of such parameters, so they are called snapshotchanging parameters. The weighting coefficients are determined for the desired field according to the intended value of these parameters. Different techniques are developed for calculating the coefficients via snapshotchanging parameters. It should be selected regarding the number and type of the parameters. In unsteady problems, the Galerkin projection method is frequently used (Bourguet et al., 2007; Zare et al., 2018). This method can also analyze steady problems with several snapshotchanging parameters (Christ & Sattelmayer, 2018; Samadiani & Joshi, 2010; Wang et al., 2012). Rolander et al. (Rolander et al., 2006) proposed Flux Matching Procedure (FMP) for when the boundary conditions are snapshotchanging parameters. In this approach, the coefficients must result in a reconstruction that has intended boundary fluxes. Rambo and Joshi (Rambo & Joshi, 2007) reported more robustness and less accuracy for the analysis of convective systems via FMP.
For problems with several snapshotchanging parameters, the most straightforward approach for the 's is the family of interpolation/extrapolation methods (Du et al., 2013; Qamar & Sanghi, 2009; Selimefendigil, 2013). It becomes a simple onedimensional curvefitting when there is only one parameter. In this family, the unknown coefficients are approximated between known coefficients obtained from the reconstruction of snapshots. The computational cost and accuracy are affected by the estimation procedure used. Wang et al. (Wang et al., 2012) compared some interpolation methods and the Galerkin projection. The Kriging (KRG) method was primarily developed for geostatistics surveys (Davis, 2002; Webster & Oliver, 2007). It is also applied in other highdimensional engineering models (Bouhlel & Martins, 2019) and thermofluid optimizations (Husain et al., 2013). Rolander suggested the Kriging method for the coefficient approximation during POD analysis (Rolander, 2005), but no one has investigated this idea yet.
In this paper, a vectorvalued POD framework is introduced for loworder, costeffective modeling of a platefin MSHE. The coefficient prediction is presented as a spatial estimation problem. A comparison is made between two methods with different degrees of complexity, i.e., the KRG and Inverse Distance Weighted (IDW) methods. To the best of our knowledge, no published article deals with such an investigation for a POD problem. The threestream MSHE is modeled numerically to collect the snapshots. Three fluid streams and the solid body of the intended platefin MSHE are analyzed conjugately. Snapshotchanging parameters are the Reynolds numbers of the flow streams. Finally, the total heat exchange, temperature distribution, and parasitic longitudinal conduction through the solid body are extracted from the reconstructions. The study steps are shown schematically in Fig. 1.
Snapshot modeling is completely according to the analysis of Shafieialamooti and Ashrafizadeh (Shafieialamooti & Ashrafizadeh, 2021). They comprise geometry and dimensions, boundary conditions, and computational domain. The threedimensional geometry and the computational domain are shown in Fig. 2(a) and (b). It is a sample design used in energy integration systems (Mehrpooya et al., 2017) and LNG plant Cold Box (Zargoushi et al., 2020). S1, S2, and S3 denote initially cold, intermediate, and hot streams. The inlet temperatures are 290 K, 330 K, and 370 K. The dimensions and thicknesses are as follows (Shafieialamooti & Ashrafizadeh, 2021): H=2.00 (mm), δ=0.10 (mm), S=1.40 (mm), and L=306.00 (mm). Flow passages are mini channels with a hydraulic diameter of 1.54 (mm) (Kandlikar & Grande, 2003).
Grid generation is also similar to that of Shafieialamooti and Ashrafizadeh (Shafieialamooti & Ashrafizadeh, 2021), but it is finer in this paper because of special postprocessing requirements. Two flow arrangements are considered for the analysis, denoted PNN and PNP, as shown in Fig. 3. In this figure, two symbols and are representatives of inward and outward directions relative to the image plane. The flow arrangement specifications PNN and PNP show the flow direction of S1, S2, and S3 streams: "P" corresponds to the flow in the positive direction along the zaxis, and "N" corresponds to the negative one. PNN and PNP arrangements are the same as P4 and P2 of Sekulic and Shah (Sekulić & Shah, 1995), respectively. They are also similar to case 3 and case 4 of Krishna et al. (Krishna et al., 2013).
Snapshots are collected by the steadystate solution of governing equations via ANSYS Fluent commercial software. Necessary grid independency tests and evaluation of numerical results are done, as explained by Shafieialamooti and Ashrafizadeh (Shafieialamooti & Ashrafizadeh, 2021). The evaluation was done during two comparisons showing a fair agreement: first with the results of Kays and London (Kays & London, 1984), and then with the temperature distribution reported by Krishna et al. (Krishna et al., 2013). The results showed good agreement according to the average difference between experimental and numerical results in comparison with similar research (Bala Sundar Rao et al., 2013; Chu et al., 2014; Yin & Ooka, 2015).
The MSHE is modeled considering the compressible flow of air with temperaturedependent thermophysical properties. Streams Reynolds numbers ( ) are selected to cover the laminar flow domain, so they are , and . The maximum number of combinations for PNN or PNP flow arrangements is , but a few snapshots are selected due to some limitations in the calculation process. Selected snapshots are all 27 combinations of , and , plus all eight combinations of and . The total number of snapshots for a flow arrangement will be .
The snapshots are calculated using a cluster of GHz CPUs with 50 GB of RAM. Convergence was attained after sec iterations for each snapshot on average.
M independent snapshot vectors are shown as .... . For a domain with N nodes, i.e., N degree of freedom, a general vector of will be as follows:
(1) 

In this paper, the vectorvalued approach is implemented, so each component of a snapshot vector is a vector of five dimensionless variables by itself:
(2) 

The variables in the are dimensionless versions of the primitive variables. For example:
(3) 

The variables , , , and are calculated similarly. They are projected into the range of 0.01.0 to eliminate the inconsistency in the dimension and scale of the variables. In this study, is defined as the affine subspace and analyzed via POD (Rowley & Marsden, 2000). vectors show the fluctuations around the mean values.
An inner product is defined concerning dimensionless variables to organize a vector space:
(4) 

is the volume corresponding to the i^{th} node. The geometric volume, , is a regular choice for . In general, may be redefined according to the field definition rather than x, y, and z.
The required set of orthogonal base vectors, 's, are extracted from the snapshots via POD analysis (Holmes et al., 1998; Sirovich, 1987). A linear combination of the bases will reconstruct the field for the desired situation, :
(5) 

In Eq. (5), Coefficient weights the k^{th} base. In the Method of Snapshots proposed by Sirovich (Sirovich, 1987), the bases are described according to the snapshot vectors:
(6) 

is the coefficient of the j^{th} snapshot for describing the k^{th} base. The best set of bases is calculated via the following eigenvalue problem achieved from the calculus of variations (Sirovich, 1987):
(7) 

is the correlation matrix, with the elements, 's, (Sirovich, 1987):
(8) 

It is a common practice to number the eigenvalues of in descending order and so for the corresponding eigenvectors. All of the eigenvectors are compiled into the matrix, so the elements of matrix are the coefficients of Eq. (6). In other words, the j^{th} component of the k^{th} eigenvector is (Sirovich, 1987).
As could be understood from Eq. (8), M orthogonal bases will be extracted for the correlation matrix , which is often less than N. On the other hand, the necessary number of bases for an acceptable reconstruction, m, is even less than M; thus, there are only a few unknown coefficients for a reconstruction. This leads to a significant reduction in computations and is the main advantage of POD analysis as a ReducedOrder Method (ROM).
The parameters which assess each base's importance are the eigenvalues, 's (Holmes et al., 1998; Sirovich, 1987). Relative Information Content (RIC) (Christ & Sattelmayer, 2018), , is defined as an index for the influence of first n bases in the reconstruction, as follows:
(9) 

In this paper, the weighting coefficients in Eq. (5), 's, are determined by interpolation/extrapolation methods, as discussed in the next section.
Eq. (6), as a base for the coefficient estimation, could be rewritten as:
(10) 

In other words, the bases reconstruct the affine subspace in the form of:
(11) 

So contains the necessary data for the estimation of the desired coefficients. is the coefficient of the i^{th} base vector, , for reconstructing the j^{th} affine subspace vector, .
The calculation process of the coefficients is done by Tecplot 360 2018 software. It is treated as a spatial estimation process according to three snapshotchanging parameters, , as coordinates. Kriging (KRG) and InverseDistance Weighted (IDW) techniques are examined, and their results are compared.
The Kriging (KRG) method utilizes the realization of stochastic processes as a model for unknown functions (Webster & Oliver, 2007). Here its universal formulation is implemented, and the coefficients are calculated as (Martin & Simpson, 2005; Webster & Oliver, 2007). , known as the drift, is the overall trend of the coefficients data. Smallscale variations are covered by adding the departure from the drift, .
A linear polynomial function is used for estimating . Approach for is based on an unbiased estimation with a minimum variance of errors (Webster & Oliver, 2007). An appropriate Spatial Correlation Function (SCF) must be assumed for the spatial distribution of covariance, . Spherical SFC, as an ideal form, is used for (Davis, 2002) here:
(12) 

The distance, , is relative to the s^{th} known coefficient:
(13) 

is the variance of the known coefficients. is the actual range, assumed to be percent of the Reynolds numbers’ full range (Davis, 2002):
(14) 

Various percentages were examined, and the value of was selected according to the best result. Regarding the chosen values of Reynolds numbers explained in section (3), the actual range is .
For the InverseDistance Weighted (IDW) method, the desired coefficient, , is calculated by the following function of the inverse of distances from the known ones (Webster & Oliver, 2007):
(15) 

The known coefficient, , is attributed to the reconstruction of the s^{th} snapshot. is the weighting exponent, and is the number of the nearest considered points. The results are achieved with the values of and .
Two sets of 35 bases are calculated via Intel Math Kernel Library (MKL) 11.2 according to the PNN and PNP flow arrangements. The time elapsed for extracting all 35 bases for each flow arrangement is about 10800 seconds, less than half of that taken for a converged snapshot solution.
Fig. 4 demonstrates the relative information content, E_{1n}. As shown, for , eight and nine bases are needed for the PNN and PNP arrangement, respectively; thus, it could be concluded that the necessary number of bases is not substantially affected by the flow arrangement.
The Streams Reynolds numbers chosen for the reconstruction are two sets. No. 1 has , and . No. 2 has equal Reynolds numbers, . In no. 1, is out of the range considered for snapshots, i.e., . This value necessitates a partial extrapolation, so the model's accuracy could be examined in comparison to no. 1 with full interpolation. The total number of reconstructions for all PNN and PNP arrangements is : PNN1, PNN2, PNP1, and PNP2. All Reconstructions using all 35 bases are accessed via (ShafieiAlamooti & Ashrafizadeh, 2020).
Two sets of parameters are provided and studied using the reconstructed thermal field:
The fields are reconstructed using the different numbers of bases, and the results are compared. Corresponding CFD results are the references for assessing the deviation. The postprocessing is done by Tecplot 360 2018 software.
The total heat exchanged by each of the streams, , for PNN1 and PNN2 cases is presented in Fig. 5. It is calculated from the reconstructions using various numbers of bases. Fig. 6 shows their exit temperatures, . The corresponding diagrams for the cases PNP1 and PNP2 are in Figs. 7 and 8. is the massweighted mean value over the exit section.
Relative deviation percentages for these two parameters are defined as:
(16) 

(17) 

and are the differences between the reconstruction values and the CFD ones. is calculated regarding the balance of heat conductions from upper and lower separating plates:
(18) 

By defining , the temperatures scale is eliminated for a better comparison:
(19) 

In which is the massweighted mean temperature in the stream exit section:
(20) 

Bar charts of and for both PNN and PNP arrangements are illustrated in Fig. 9. Also, Table 1 shows their average values over all three streams.
As expected and also formerly reported (Du et al., 2013), there is a negligible deviation in the cases of PNN2 and PNP2 due to the full interpolation. In the cases of PNN1 and PNP1, the deviations are significant for some of the results because of the partial extrapolation.
Figs. 58 depict that increasing the number of bases up to 15 may decrease the deviation. This number of bases is more than that predicted in Fig. 4 according to . The postprocessing integrations for and engage multiple variables according to Eq. (18) and (20). This leads to the integration of the deviation for several variables (Gerald & Wheatley, 2004). The increased number of bases is to compensate for it. In addition, the deviation of temperature causes a secondary deviation due to the temperaturedependent variables in Eq. (18) and (20), i.e., , , and .
Fig. 9 shows a higher deviation for in comparison with . From Table 1, the average of is 6.2 for the IDW and 7.2 for the KRG. This observation could be due to the larger area of integration in Eq. (18) relative to Eq. (20), which allows for a higher deviation for .
Figs. 58 reveal that the KRG method often has better results in comparison with the IDW method. KRG method outperforms the IDW method in the full interpolations, i.e., PNN2 and PNP2 cases. The mean ratio of deviations for these interpolated cases, , is 9.3, and the deviations for both methods are negligible in some cases.
The IDW method has less deviation in four resulting parameters out of the total 24 ones, i.e., 17% of them. These results are and for the PNN1 and PNP1 cases. In these two cases, the deviation of is again lower for the KRG. Both methods have nonnegligible deviations for in PNP1, about 96% with KRG and 160% with IDW; thus, about 4% of all estimations are unacceptable. It could be recognized that for PNP1 is the smallest value of between three streams in all cases. Recalling that it is the denominator of Eq. (16), the severe rise of could be described. The figures clarify a more substantial contribution of intermediate stream in the whole heat exchange process for the PNN cases relative to the PNP ones. It could be concluded that the POD has a weaker capability to estimate a nondominant phenomenon.
According to Fig. 9, deviations could be different for the cases with similar sets of but different flow directions; thus, the flow arrangement is a determinant ofthe usability of POD analysis.
Longitudinal temperature distributions are illustrated in Figs. 1013 for the cases PNN1, PNN2, PNP1, and PNP2, respectively. They are shown as vs. . The is calculated according to Eq. (19). The results are approximated by 1, 3, 5, 7, 9, and all 35 bases.
For the PNN2 and PNP2, the POD curves are closely approaching the CFD curves by a few bases. For the PNN1 and PNP1, the trends are predicted reasonably for the cold and hot streams, but the estimations are not so satisfactory for the intermediate stream. There are three extrema in the temperature distribution curves, all of which are for the intermediate stream. Only for the PNP1, the position of the extremum point has a nonnegligible displacement relative to the CFD prediction. This shifting is observed for both KRG and IDW methods. In Fig. 10, for the curve of the cold stream for PNN1, the maximum deviation occurred in different positions for KRG and IDW: there is more deviation in the middle parts for KRG, while for the IDW the deviation is considerable in the exit section, i.e., .
Approximations for temperature contours are shown in Figs. 14 and 15. The figures clarify a more substantial contribution of the intermediate stream in the whole heat exchange process for the PNN flow arrangement relative to the PNP. The difference between the two methods is distinguishable only for the PNP1. Referring to the deviation values, it could be concluded that the POD has a weak capability for the estimation of a nondominant phenomenon.
Longitudinal conductive heat transfer through the solid body is one of the parasitic effects of the thermal process. Here it is studied in terms of vs. , and shown in Fig. 16. is the value of longitudinal conduction, and is the maximum heat exchange between streams.
In the PNN2 and PNP2 cases, KRG and IDW predictions are satisfactory using 2 and 4 bases, respectively. In PNN1 and PNP1, both methods cover the trends qualitatively. KRG generally has better estimations, but there are exceptions: For the PNN1, there are two minimum points near two ends and IDW has a better prediction for the minimum value near .
The temperature distribution in the solid body governs the longitudinal conduction phenomenon. It is mainly the function of two parameters: the temperatures of three streams as shown in Figs. 1013, and their heat capacity rates (Shafieialamooti & Ashrafizadeh, 2021). It is seen from Figs. 1013 that there are higher deviations for temperatures in the middle of the stream passages leading to a more significant deviation of there. On the other hand, a severe temperature deviation for the intermediate stream does not corrupt the result for because of the contribution of the other two streams in the prediction of .
A threestream platefin heat exchanger is analyzed by proper orthogonal decomposition. Snapshots are collected from a conjugate numerical simulation of fluids and solid bodies. Bases are extracted using a vectorvalued approach with dimensionless variables. Coefficient prediction is treated as a spatial estimation procedure, and two methods of Kriging (KRG) and InverseDistance Weighted (IDW) are examined.
Two flow arrangements are considered in which two combinations of the Reynolds numbers are reconstructed. In the first one, all three Reynolds numbers, Re^{cold}, Re^{int}, and Re^{hot}, are in the snapshots’ domain, so the coefficients are fully interpolated. In the second one, Re^{int} is out of the domain; thus, it requires partial extrapolation. The results are compared with the CFD predictions. The main conclusions are as follows:
Fig. 1. Schematic diagram of the study steps.
Fig. 4. The relative information content, E_(1n).
Table 1. Average values of the relative deviation percentages over all three streams
Arrangement Case no. ϵ_Q^ave ϵ_T^ave
InverseDistance Kriging InverseDistance Kriging
PNN 1 15.43 12.98 7.32 5.78
2 2.11 0.31 0.74 0.13
PNP 1 58.17 37.69 4.03 1.24
2 2.46 0.42 0.56 0.03