A POD Approach for the Analysis of a Multi-Stream Heat Exchanger of a Cold Box: Comparison of two different Coefficient Approximations

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 (Shafiei-Alamooti & Ashrafizadeh, 2020).

Article Highlights

• Variable-property conjugate modeling is used to collect snapshots.
• Base coefficients are approximated as spatial estimations.
• The Kriging coefficient approximation has better results in 83% of cases.
• The necessary number of bases depends on the required post-processing.

Acknowledgments

Funding: This work was supported by the National Iranian Gas Company, NIGC [grant number 950406].

1. Introduction

Multi-Stream 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 two-stream ones; thus, there is a high demand for cost-effective 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. Phase-change 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 multi-stream 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 large-scale 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). Shafiei-alamooti and Ashrafizadeh (Shafiei-alamooti & Ashrafizadeh, 2021) did a phenomenological analysis on a plate-fin 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 cost-effective techniques. Soft computing methods like Neural Networks and Genetic algorithms are frequently used in conventional two-stream 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 plate-fin heat exchanger. Yin and Ooka (Yin & Ooka, 2015) simulated a passage of a plate-fin 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 thermo-fluid systems, such as data center server cooling systems (Ghosh & Joshi, 2014), automotive HVAC systems (Christ & Sattelmayer, 2018), thermo-acoustic heat engine (Selimefendigil & Öztop, 2014), and Phase-Change 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 air-cooled 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 POD-aided 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):

1. Collection of the known data, numerically or experimentally, to be used as snapshots.
2. Extraction of the orthogonal bases by an eigenanalysis of snapshot data.
3. Searching appropriate weighting coefficients of the bases, 's, for the intended unknown field.

For a thermo-fluid problem, POD analysis is strongly affected by the selection of the configuration of orthogonal bases: the scalar-valued approach (Brenner et al., 2009; Rowley et al., 2001) or the vector-valued 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 snapshot-changing 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 snapshot-changing 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 snapshot-changing 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 snapshot-changing 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 snapshot-changing 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 one-dimensional curve-fitting 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 high-dimensional engineering models (Bouhlel & Martins, 2019) and thermo-fluid 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 vector-valued POD framework is introduced for low-order, cost-effective modeling of a plate-fin 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 three-stream MSHE is modeled numerically to collect the snapshots. Three fluid streams and the solid body of the intended plate-fin MSHE are analyzed conjugately. Snapshot-changing 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.

2. Thermo-Hydraulic Model of Snapshots

Snapshot modeling is completely according to the analysis of Shafiei-alamooti and Ashrafizadeh (Shafiei-alamooti & Ashrafizadeh, 2021). They comprise geometry and dimensions, boundary conditions, and computational domain. The three-dimensional 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 (Shafiei-alamooti & 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 Shafiei-alamooti and Ashrafizadeh (Shafiei-alamooti & Ashrafizadeh, 2021), but it is finer in this paper because of special post-processing 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 z-axis, 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 steady-state solution of governing equations via ANSYS Fluent commercial software. Necessary grid independency tests and evaluation of numerical results are done, as explained by Shafiei-alamooti and Ashrafizadeh (Shafiei-alamooti & 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 temperature-dependent thermo-physical 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.

3. POD-Based Reduced - Order Reconstruction Method
   • Orthogonal Bases Extraction

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:

In this paper, the vector-valued approach is implemented, so each component  of a snapshot vector  is a vector of five dimensionless variables by itself: 

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

The variables , , , and  are calculated similarly. They are projected into the range of 0.0-1.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:

 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, :

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:

 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):

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

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 Reduced-Order 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:

In this paper, the weighting coefficients in Eq. (5), 's, are determined by interpolation/extrapolation methods, as discussed in the next section.

• Approximation of Weighting Coefficients

Eq. (6), as a base for the coefficient estimation, could be rewritten as:

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

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 snapshot-changing parameters, , as coordinates. Kriging (KRG) and Inverse-Distance 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. Small-scale 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:

The distance, , is relative to the s^th known coefficient:

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

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 Inverse-Distance Weighted (IDW) method, the desired coefficient, , is calculated by the following function of the inverse of distances from the known ones (Webster & Oliver, 2007):

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 .

4. Results and Discussion

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_1-n. 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 : PNN-1, PNN-2, PNP-1, and PNP-2. All Reconstructions using all 35 bases are accessed via (Shafiei-Alamooti & Ashrafizadeh, 2020).

Two sets of parameters are provided and studied using the reconstructed thermal field:

• Major quantities show the general performance of the thermal process, including streams' heat exchange and exit temperature.
• Minor quantities; consist of the distribution of the temperatures along the flow paths and the longitudinal heat conduction in the solid body of MSHE.

 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 post-processing is done by Tecplot 360 2018 software.

• Major Quantities: Heat Exchange and Exit Temperature

The total heat exchanged by each of the streams, , for PNN-1 and PNN-2 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 PNP-1 and PNP-2 are in Figs. 7 and 8.  is the mass-weighted mean value over the exit section.

Relative deviation percentages for these two parameters are defined as:

 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:

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

 In which  is the mass-weighted mean temperature in the stream exit section:

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 PNN-2 and PNP-2 due to the full interpolation. In the cases of PNN-1 and PNP-1, the deviations are significant for some of the results because of the partial extrapolation.

Figs. 5-8 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 post-processing 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 temperature-dependent 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. 5-8 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., PNN-2 and PNP-2 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 PNN-1 and PNP-1 cases. In these two cases, the deviation of  is again lower for the KRG. Both methods have non-negligible deviations for  in PNP-1, about 96% with KRG and 160% with IDW; thus, about 4% of all estimations are unacceptable. It could be recognized that  for PNP-1 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 non-dominant 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.

• Minor Quantities
  • Temperature Distribution

Longitudinal temperature distributions are illustrated in Figs. 10-13 for the cases PNN-1, PNN-2, PNP-1, and PNP-2, 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 PNN-2 and PNP-2, the POD curves are closely approaching the CFD curves by a few bases. For the PNN-1 and PNP-1, 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 PNP-1, the position of the extremum point has a non-negligible 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 PNN-1, 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 PNP-1. Referring to the deviation values, it could be concluded that the POD has a weak capability for the estimation of a non-dominant phenomenon.

• Longitudinal Conduction

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 PNN-2 and PNP-2 cases, KRG and IDW predictions are satisfactory using 2 and 4 bases, respectively. In PNN-1 and PNP-1, both methods cover the trends qualitatively. KRG generally has better estimations, but there are exceptions: For the PNN-1, 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. 10-13, and their heat capacity rates (Shafiei-alamooti & Ashrafizadeh, 2021). It is seen from Figs. 10-13 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 .

5. Conclusion

A three-stream plate-fin 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 vector-valued approach with dimensionless variables. Coefficient prediction is treated as a spatial estimation procedure, and two methods of Kriging (KRG) and Inverse-Distance 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:

1. For the major quantities, i.e., total heat exchange and exit temperature, the KRG outperforms the IDW for 87% of the results. For the cases with full interpolation, accurate estimations are achieved by only a few bases. In these cases, the average deviation of KRG is 9.3 times lower than the IDW. The estimations' accuracy falls when a partial extrapolation is needed. In these cases, IDW has less deviation in 17% of estimated parameters.
2. The necessary number of bases for the minimum deviation depends on the calculation process of the intended parameter and may be higher than those predicted by the eigenvalues. Interaction between the POD analysis and the property variation may lead to deviation growth. Also, the flow arrangement could affect the deviation, so it may be decisive for the applicability of unified POD analysis.
3. Longitudinal heat conduction is governed by the temperatures of all three streams, so its deviation rises when there are poor approximations for the temperatures of streams. On the other hand, a deviation from the extrapolation for one stream may not be so corruptive because of the effect of two other streams.
4. Comparison of the total heat exchanges with their deviations reveals that the POD has a weak capability for estimating non-dominant phenomena and infinitesimal values.
5. KRG method utilizes a more extensive set of known data for the estimation, so it has almost better results, as well as a higher computational cost. However, the results may be improved for both KRG and IDW by further detailed investigations about the parameters involved.
6. The optimized design of cold boxes could involve several trial-and-errors, so the suggested method could facilitate and advance it by reducing the computational cost and time. The method's performance strongly depends on the accuracy of the snapshot data and the domain of snapshot collection. Finally, a sufficient number of snapshots could be collected by detailed numerical modeling for a particular industrial arrangement consisting of phase-change phenomena. Also, possible experimental research could support both snapshot collection and coefficient calculation.

 
Fig. 1. Schematic diagram of the study steps.
  
Fig. 4. The relative information content, E_(1-n).

  
Table 1. Average values of the relative deviation percentages over all three streams
Arrangement    Case no.    ϵ_Q^ave    ϵ_T^ave
        Inverse-Distance    Kriging    Inverse-Distance    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