Analysis of the Mathematical Models for Identifying the Thickness of the Fouling Layer in Natural Gas Coolers

Abstract

This article presents an analysis of three different approaches to the identification of the thickness of the fouling layer inside the pipes of natural gas (NG) coolers. At present, there is no existing simple analytical procedure for the identification of the fouling layer thickness. The authors of this article describe in detail the balance method, which required the use of a large number of physical parameters, changes in their sizes depending on the output temperature of the gas, the temperature of the cooling air, the air quantity, as well as the physical properties of both media. The computational model was robust, and its disadvantage was the iterative computation. The second analysed method was a dimensional analysis. It was applied using the Buckingham’s theorem to express the individual similarity criteria. In this method, 10 simplexes and two complexes were created. The fouling layer thickness, expressed using a derived criterial equation, exhibited real results. The third analysed method was based on analysing selected physical parameters with the use of a multiple regression analysis in MinitabX 18 software. The analysis showed that the fouling layer thickness depended on fewer parameters than the number of parameters assumed in the dimensional analysis or the balance method. The standard deviation that was identified in the multiple linear regression for a double crossflow cooler was 0.0667 and the value of reliability (the coefficient of determination of the multiple linear regression) R² was 0.9985.

1. Introduction

Natural gas coolers are pressure vessels with many thousands of sealed holes. In the coolers, cooled natural gas flows through tubes with external fins. The maximum projected gas pressure in coolers of various designs is 7.85 MPa. During the cooler operation, a fouling layer is formed on the inside wall of the tubes, and it reduces the cross-sectional area of the tubes. As the fouling layer thickness increases, several negative effects occur: a higher speed of the flowing gas, higher thermal resistance of the tube wall, changes in the operating conditions, etc.

The motivation for conducting the study presented in this article was to find a simple procedure that would facilitate the identification of the thickness of a fouling layer formed on a heat-transfer surface from any material whatsoever. Such materials may include, for example, a limescale layer inside the pipes of a heat exchanger, sediments on the walls of a plate heat exchanger in the pasteurisation of soy milk, and so forth. The currently available literature only provides general equations describing the cooling process for technical devices that require extensive iteration procedures. The authors of the present article have developed a simpler procedure that facilitates the identification of a real thickness of the fouling layer formed on the internal heat-transfer surfaces in the pipes of a cooler. Therefore, this study includes a rather detailed analysis of mathematical equations. The proposed method brings relevant results that have been validated in experiments.

Natural gas coolers are characterised by cooling surfaces of various sizes and shapes. The numbers of tubes and their arrangements inside the cooler block, as well as the methods of transporting gas through the tubes, may vary. Over a prolonged operation of a cooler, its cooling capacity decreases due to sediments accumulated on the internal and external surfaces of the tubes. Keeping the external heat-transfer surface clean is relatively easy; however, cleaning the internal surfaces is a complicated process. The cooler must be dismantled, all its tubes must be cleaned with the use of special technology, and then it must be re-assembled. As the cleaning technology is complicated and very long, the cooler must be removed from service for a long period of time; this causes a considerable financial loss to the gas transporter.

The currently available analytical procedures for the identification of a cooling capacity of any equipment in which heat is transferred between the cooled medium and the cooling medium are rather complicated. In general, the key challenge is to exactly identify the heat transfer coefficient for both the cooled medium and the cooling medium, especially in cases where the internal or external heat-transfer surfaces are fouled and the external surface has a complex shape. Figure 1 shows the real condition of a natural gas cooler fouled with impurities of an unknown composition. The authors of paper [1] identified the heat-transfer coefficient of the fouling layer in an experiment. Expressing a cooling capacity is also difficult due to the fact that the temperatures of the cooled medium and of the cooling medium change on the heat-transfer surface along the entire tube length [2,3,4]. The basics of the process of cooling gas media by applying ventilation technology may be found in the literature [5,6]. The analysed cooler was designed with the tubes arranged in several rows on top of each other. In each row, the temperatures of the gas and of the cooling air are different. As the temperatures change, the heat transfer coefficients for the internal and external sides also change, and so does the gas density and viscosity, as well as the thermal conductivity coefficient.

The authors were searching for a procedure that would describe in the simplest possible way, but with sufficient accuracy, the process of formation of the fouling layer in natural gas coolers that are used in various compressor stations and have various shapes of the heat-transfer surface. In this article, two different methods for the identification of the fouling layer thickness were analysed: the balance method, and the procedure using dimensional analysis. When applying dimensional analysis, it is advisable to confront the obtained results with the results obtained by a physical experiment [7] or numerical modelling [8,9], or a complicated analytical procedure [10,11,12]. The method that is described in this article is based on the basics of the probability theory and modelling presented in papers [13,14,15,16,17,18]. Dimensional analysis was also applied in papers [19,20].

The authors of this article have published their own studies on the application of dimensional analysis in various fields in relevant journal studies and publications. A summary of the modelling methods applied with the use of dimensional analysis is presented in publication [21].

For the purpose of proving the applicability of the created mathematical model, the task presented in this article was carried out with a double crossflow cooler. Such coolers are used in compressor stations, for example in Veľké Kapušany, Veľké Zlievce, and elsewhere in Slovakia, as well as in other countries.

Each of the aforesaid methods is described in detail in Section 2.

2. Mathematical Models for Expressing the Fouling Layer Thickness

Each of the methods described below has its positive and negative sides. It is important that the used method is as simple as possible in expressing the fouling layer thickness and provides real results compared to the complicated procedures. In this study, the experimental verification of results was not acceptable for the reasons mentioned in the introductory and final sections of this article. Therefore, the comparison will involve the results obtained by applying the individual methods described below.

2.1. Balance Method

The balance method is based on a comparison of the supplied and removed heat output using Equation (1). The method applicable to a crossflow cooler (Figure 2a) is described in detail in paper [1], while the results related to a double crossflow cooler (Figure 2b) are presented in this article.

The equation for calculating the cooling capacity P is as follows:

$$P=P_1·{\eta }_{\mathrm{l}\mathrm{o}\mathrm{s}}=Q_{m,\mathrm{N}\mathrm{G}}·\left(i_{1,\mathrm{N}\mathrm{G}}-i_{2,\mathrm{N}\mathrm{G}}\right)·{\eta }_{\mathrm{l}\mathrm{o}\mathrm{s}}=Q_{m,\mathrm{a}}·\left(i_{2,\mathrm{a}}-i_{1,\mathrm{a}}\right)$$

(1)

wherein P₁ is the heat output withdrawn from natural gas (W); η_los is the coefficient respecting the heat loss into the cooler’s structure (1); Q_m,NG or Q_m,a is the mass flow rate of gas or air flowing through the cooler (kg∙s⁻¹); i_1,NG, i_2,NG is the specific enthalpy of gas at the entry into and exit from the cooler (J∙kg⁻¹); and i_1,a, i_2,a is the specific enthalpy of air at the entry into and exit from the cooler (J∙kg⁻¹).

Assuming that, for example, 3% of the output P₁ transfers into the cooler’s structure, then coefficient η_los equals 0.97 and the following equation is applied:

$$P={0.97·P}_1$$

(2)

Concurrently with Equation (1), the following equation is applied:

$$P=k·S_{\Sigma }·\overline{∆t}$$

(3)

wherein k is the overall heat transfer coefficient (W·m⁻²·K⁻¹); $\overline{∆t}$ is the mean temperature gradient (K); and S_Σ is the total heat-transfer surface area of the cooler (m²).

The total surface area of the cooler S_Σ in Equation (3) is determined by the product of the external surface area of a single tube S_e (m²) and the number of all tubes in the cooler n_p:

$$S_{\Sigma }=n_{\mathrm{p}}·S_{\mathrm{e}}$$

(4)

The surface area of a single tube is calculated using the following equation:

$$S_{\mathrm{e}}=r·\pi \left(\frac{1}{2}·(d_{\mathrm{r}}^2-d_2^2)+d_2·b\right)$$

(5)

wherein r is the number of fins on a single tube (1); d_r is the fin diameter (m); d₂ is the external diameter of the tube (m); and b is the fin spacing (m).

The mean temperature gradient $\overline{∆t}$ (K) in Equation (3) may be identified using the equations presented in detail in paper [1]. The analysed cooler used a double crossflow with two flow runs n = 2. The correction coefficient ψ that is necessary for the identification of the mean temperature gradient $\overline{∆t}$ is therefore calculated using the following equation:

$$\psi =\frac{\mathrm{l}\mathrm{n}A}{n·(1-\mathrm{R})·\mathrm{l}\mathrm{n}\left(1+\frac{1}{R}·\mathrm{l}\mathrm{n}\frac{R-1}{R·A^{\frac{1}{n}}-1}\right)}$$

(6)

For parameter A, the Equation (6) determines the following:

$$A=\frac{1-\mathrm{P}}{1-R·\mathrm{P}}$$

(7)

The P and R criteria are identified using the following equations:

$$P=\frac{K_{\mathrm{N}\mathrm{G}}}{{\mathrm{K}}_{\mathrm{N}\mathrm{G}}+K_{\mathrm{a}}} ; R=\frac{K_{\mathrm{N}\mathrm{G}}·{\eta }_{\mathrm{l}\mathrm{o}\mathrm{s}}}{K_{\mathrm{a}}}$$

(8)

wherein K_NG is the heat capacity rate of gas; K_a is the heat capacity rate of air (W·K^-1).

The heat capacity rate expresses the output that is required to heat the mass of the flowing medium by 1 K and is calculated using the following equation:

$$K=Q_m·c_p$$

(9)

wherein c_p is the specific heat capacity of the medium (J∙kg⁻¹∙K⁻¹).

After substituting the total heat-transfer surface area and the mean temperature gradient into Equation (3), the overall heat transfer coefficient k may be calculated using Equation (10):

$$k=\frac{P}{S_{\Sigma }·\overline{∆t}}$$

(10)

Coefficient k may also be calculated based on a detailed description of the process of heat transport from warm natural gas through the tube wall into the surrounding air. If there is a fouling layer inside the tube, the solution must also take into account the heat transfer through that fouling layer.

The authors of this article assumed that the cooler had a fouling layer inside and a clean external heat-transfer surface, and they derived the following equation for calculating the overall heat transfer coefficient:

$$k=\frac{1}{\left(\frac{1}{{\alpha }_1}+\frac{h_{\mathrm{f}}}{{\lambda }_{\mathrm{f}}}\right)·\frac{S_{\mathrm{e}}}{\pi ·\left(d_1-2·h_{\mathrm{f}}\right)·L}+\frac{s_{\mathrm{p}}}{{\lambda }_{\mathrm{p}}}·\frac{S_{\mathrm{e}}}{S_1}+\frac{1}{{\alpha }_2}}$$

(11)

wherein α₁ is the coefficient of heat transfer inside the tube (W∙m⁻²∙K⁻¹); h_f is the thickness of the fouling layer in the tube (m); λ_f is the thermal conductivity coefficient of the fouling layer (W∙m⁻¹∙K⁻¹); λ_p is the thermal conductivity coefficient of the tube (W∙m⁻¹∙K⁻¹); S₁ is the internal surface area of the tube (m²); d₁ is the internal diameter of the tube (m); s_p is the thickness of the tube wall; L is the tube length (m); and α₂ is the heat transfer coefficient for the heat transfer from the finned tube into the air (W∙m⁻²∙K⁻¹).

By comparing Equations (11) and (10), the following equation for calculating the fouling layer thickness was obtained:

$$h_{\mathrm{f}}=\frac{D·\pi ·L·d_1·\frac{S_{\mathrm{e}}}{{\alpha }_1}}{\frac{S_{\mathrm{e}}}{{\lambda }_{\mathrm{f}}}+2·D·\pi ·L}$$

(12)

Parameter D in Equation (12) was as follows:

$$D=\frac{S_{\Sigma }·\overline{∆t}}{P}-\frac{s_{\mathrm{p}}}{{\lambda }_{\mathrm{p}}}·\frac{S_{\mathrm{e}}}{S_1}-\frac{1}{{\alpha }_2}$$

(13)

The key challenge associated with calculating the fouling layer thickness using Equation (12) was that some of the parameters in Equation (12) are a function of h_f. As h_f increases, the thermal resistance on the tube wall also increases, and this leads to a decrease in the cooling capacity P. The fouling layer reduces the internal diameter of the tube and, consequently, the natural gas flow rate increases as well as the value of coefficient α₁. The heat transfer coefficient α₁ is calculated using a criterial equation that contains, inter alia, the Reynolds number Re. A higher α₁ value leads to a lower value of thermal resistance on the interface between the natural gas and the internal wall of the tube. As a result, the cooling capacity increases. Moreover, changes in the fouling layer thickness cause changes in the mean temperature gradient $\overline{∆t}$. Based on the aforementioned facts, it may be stated that the fouling layer thickness cannot be calculated directly and requires numerous iterations.

The basic geometric parameters of the analysed cooler are listed in

Table 1. Geometric parameters of the analysed cooler.

Parameter
d1	d2	dr	sr	B	L	s1	np
m	m	m	m	M	m	m	1
0.025	0.033	0.057	0.0007	0.0026	9	0.064	2268

Parameter s_r is the fin thickness and s₁ is the tube spacing.

.

The efficiency of the gas cooling process is significantly affected by the thermal conductivity of the fouling layer. The authors of this article identified the thermal conductivity coefficient of the fouling layer in an experiment. The coefficient value was 0.746 W·m⁻¹·K⁻¹, approximately two orders lower than the value for the tube material [1].

The key settings for the iterations were selected based on the project values of the cooler’s parameters. The cooler was designed to cool natural gas at 242.3 kg·s⁻¹ and a pressure of 7.45 MPa from a temperature of 75 °C to 50 °C. The required quantity of the fan air with a temperature of 28 °C was 565.23 kg·s⁻¹.

In the project, it was assumed that the heat-transfer surfaces were clean; therefore, h_j = 0 was substituted into Equation (11). A similar procedure was applied to the iterations that were used for the temperatures of the cooling air which were different from the projected temperature. The ambient temperature ranged from 28 to 0 °C. The mass flow rate, the input temperature, and the gas pressure remained constant during the process. As the ambient temperature decreased, the amount of the air decreased, as described by the following equation:

$$Q_{m,\mathrm{a}}=4.8349397·{10}^{-3}·T_{\mathrm{a}}^3-3.9123122·T_{\mathrm{a}}^2+1060.1414·T_{\mathrm{a}}-95948.998$$

(14)

Selected results of iterations are presented in Figure 3. For certain values of the ambient temperature, the curve of the fouling layer thickness corresponds to the relevant gas cooling degree.

The indicated calculation method is rather complicated. Moreover, it requires obtaining extensive input data before the solving process is commenced—not only geometric data, but also the values of physical properties of natural gas and air. Particularly in the case of natural gas, with a pressure ranging from approximately 6 to 7.5 MPa, complex functions that depend on a pressure and a temperature must be applied in order to describe its physical properties.

Due to the complexity of the iteration procedure, the authors attempted to identify the fouling layer thickness through dimensional analysis based on the Buckingham’s theorem.

2.2. Buckingham’s Method

In the first approximation, all the parameters that were selected for the purpose of deriving a criterial equation were identical to those used in the balance method. In addition to the parameters listed above, these also included density, the thermal conductivity coefficient and kinematic viscosity of the flowing media, as well as the thermal conductivity coefficients of the tube and fin materials. A total of 18 dimensionless criteria were created and the criterial equation was subjected to an analysis, which revealed the effects of the individual criteria on the resulting h_f value. The parameters with apparently minor effects were then withdrawn from the model, and a new, simpler model was created. In the new model, some of the geometric dimensions (e.g., the external diameter of the tube, the fin diameter and thickness, and the fin spacing) were replaced with the surface areas calculated on the basis of those parameters. The adjusted list of physical parameters (16 in total) is presented in

Table 2. Input parameters for the mathematical model.

Symbol	Physical Parameter	Unit
KNG	Heat capacity rate of gas	kg·m2·s−3·K−1
Ka	Heat capacity rate of air	kg·m2·s−3·K−1
T1,NG	Gas temperature at the inlet into the cooler	K
T2,NG	Gas temperature at the outlet from the cooler	K
T1,a	Air temperature at the inlet into the cooler	K
S1,Σ,c	Net cross-sectional area of all tubes for gas	m2
d1	Internal diameter of the cooler tube	m
L	Tube length	m
hf	Fouling layer thickness	m
SΣ	Total heat-transfer surface area of the cooler	m2
Se	Surface area of a single finned tube in the cooler	m2
Sa,Σ	Free surface area of the cooler for air passage	m2
ρNG	Gas density	kg·m−3
ρa	Air density	kg·m−3
λNG	Thermal conductivity of gas	kg·m·s−3·K−1
λf	Thermal conductivity of the fouling layer	kg·m·s−3·K−1

.

In the creation of the model law, all the dimensions of the selected physical parameters were transformed into seven SI base units (kg; m; s; K, A; mol; cd). The relevant parameters in the aforesaid model for the identification of the fouling layer thickness included only four basic dimensions—kg, m, s, and K.

A criterial equation is normally created by replacing the selected dimensional parameters φ₁ through φ_n with similarity criteria π₁ through π_m, while the functional correlations between the individual criteria are identified experimentally or through numerical or analytical calculations. The resulting criterial equation then applies to the entire group of similar processes.

When a dimensional analysis is used to describe a process, the number of obtained criteria π is always lower than the number of relevant parameters n on which the process depends. The basic equation that expresses the corelations between n relevant parameters φ₁, φ₂ … φ_i … φ_n of various dimensions is as follows [21]:

$$f\left({\phi }_1,{\phi }_2\dots {\phi }_i\dots {\phi }_n\right)=0$$

(15)

Based on the defining equation, each of the φ parameters may be expressed through a specific dimensional equation. It is the product of the base unit symbols with the respective exponents. For the four base units of the selected physical parameters (kg, m, s, K), the defining equation is as follows:

$$\left[\phi \right]={\mathrm{m}}^{x_1}·{\mathrm{k}\mathrm{g}}^{x_2}·{\mathrm{s}}^{x_3}·{\mathrm{K}}^{x_4}$$

(16)

In Equation (16), dimensional exponents x₁ through x₄ are rational numbers that were identified as described below.

Equation (15) is dimensionally homogenous; therefore, φ_i variables cannot be used separately, but only in the form of products:

$$\pi =\prod _{i=1}^n{\phi }_i^{x_i}$$

(17)

wherein

• π—is the dimensionless variable (similarity criterion) (1);
• x_i—is the exponent (rational number);
• φ_i—is the physical parameter with a respective dimension.

Pursuant to Equation (15) and considering the physical parameters listed in Table 2, the following correlation must apply:

$$f(K_{\mathrm{N}\mathrm{G}},K_{\mathrm{a}},T_{1,\mathrm{N}\mathrm{G}},T_{2,\mathrm{N}\mathrm{G}},T_{1,\mathrm{a}},S_{1,\Sigma ,\mathrm{c}},d_1{,L,h}_{\mathrm{f}},S_{\Sigma },S_{\mathrm{e}},S_{\mathrm{a},\Sigma },{\rho }_{\mathrm{N}\mathrm{G}},{\rho }_{\mathrm{a}},{\lambda }_{\mathrm{N}\mathrm{G}},{\lambda }_{\mathrm{f}})=0$$

(18)

In general, for a certain phenomenon that is described by n relevant parameters, it is possible to create l similarity criteria. The number of searched criteria π is identified using the following equation:

$$l=n-h$$

(19)

wherein h is the je rank of the dimension matrix.

Pursuant to Equation (17), the following applies:

$$\pi =K_{\mathrm{N}\mathrm{G}}^{x_1}·K_{\mathrm{a}}^{x_2}·T_{1,\mathrm{N}\mathrm{G}}^{x_3}·T_{2,\mathrm{N}\mathrm{G}}^{x_4}·T_{2,\mathrm{a}}^{x_5}·S_{1,\Sigma ,\mathrm{c}}^{x_6}·d_1^{x_7}·L^{x_8}·h_{\mathrm{f}}^{x_9}·S_{\Sigma }^{x_{10}}·S_{\mathrm{e}}^{x_{11}}·S_{\mathrm{a},\Sigma }^{x_{12}}·{\rho }_{\mathrm{N}\mathrm{G}}^{x_{13}}·{\rho }_{\mathrm{a}}^{x_{14}}{·\lambda }_{\mathrm{N}\mathrm{G}}^{x_{15}}{·\lambda }_{\mathrm{f}}^{x_{16}}=0$$

(20)

The dimensional form of Equation (20) is as follows:

$$\begin{array}{cc}\hfill 1=(\mathrm{k}\mathrm{g}·{\mathrm{m}}^2·{\mathrm{s}}^{-3}& {·{\mathrm{K}}^{-1})}^{x_1}·({\mathrm{k}\mathrm{g}·{\mathrm{m}}^2·{\mathrm{s}}^{-3}·{\mathrm{K}}^{-1})}^{x_2}·{\mathrm{K}}^{x_3}·{\mathrm{K}}^{x_4}·{\mathrm{K}}^{x_5}·{\left({\mathrm{m}}^2\right)}^{x_6}·(\mathrm{m}{)}^{x_7}·(\mathrm{m}{)}^{x_8}·(\mathrm{m}{)}^{x_9}·{(\mathrm{m}}^2{)}^{x_{10}}·\hfill \\ & ·{\left({\mathrm{m}}^2\right)}^{x_{11}}·{\left({\mathrm{m}}^2\right)}^{x_{12}}·{(\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-3})}^{x_{13}}·{(\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-3})}^{x_{14}}·{(\mathrm{k}\mathrm{g}·\mathrm{m}·{\mathrm{s}}^{-3}·{\mathrm{K}}^{-1})}^{x_{15}}·\hfill \\ & ·{(\mathrm{k}\mathrm{g}·{\mathrm{m}}^2·{\mathrm{s}}^{-3}·{\mathrm{K}}^{-1})}^{x_{16}}\hfill \end{array}$$

(21)

Since the left side of Equation (21) equals one, the sum of the dimensional exponents in every basic parameter must equal zero. Therefore, the individual dimensions of the physical parameters (kg, m, s, K) are subject to the set of Equations (22)–(25).

For the “kilogram” unit:

$$x_1{+x}_2{+x}_{13}{+x}_{14}{+x}_{15}{+x}_{16}=0$$

(22)

For the “meter” unit:

$${2x}_1+2x_2+{2x}_6+x_7+x_8+x_9+{2x}_{10}+{2x}_{11}+{2x}_{12}-3x_{13}-{3x}_{14}+x_{15}+{2x}_{16}=0$$

(23)

For the “second” unit:

$${-3x}_1{-3x}_2{-3x}_{15}{-3x}_{16}=0$$

(24)

For the “Kelvin” unit:

$${-x}_1-x_2+x_3+x_4+x_5-x_{15}-x_{16}=0$$

(25)

The rank of the matrix for the set of Equations (22)–(25) equals four. With a total of 16 physical parameters, the total number of criteria is l = 12, as calculated using Equation (19).

Parameters with identical dimensions are expressed as independent criteria, also referred to as simplexes. The number of simplexes l_s equals the difference between the total number of relevant parameters n and the number of relevant parameters with different dimensions n_k; therefore:

$$l_{\mathrm{s}}=n-n_{\mathrm{k}}$$

(26)

Table 1 indicates that the number of parameters with different dimensions is six. This means that a task with 16 relevant parameters may be assigned 10 simplex criteria.

The simplex that was created on the basis of the heat capacity rate values, calculated for gas and air using Equation (9), is as follows:

$${\pi }_1=\frac{K_{\mathrm{N}\mathrm{G}}}{K_{\mathrm{a}}}$$

(27)

The temperature simplex criteria are as follows:

$${\pi }_2=\frac{T_{1,\mathrm{N}\mathrm{G}}}{T_{1,\mathrm{N}\mathrm{G}}-T_{2,\mathrm{N}\mathrm{G}}}=\frac{T_{1,\mathrm{N}\mathrm{G}}}{{\mathsf{\Delta }T}_{\mathrm{N}\mathrm{G}}}$$

(28)

$${\pi }_3=\frac{T_{2,\mathrm{N}\mathrm{G}}}{T_{1,\mathrm{a}}}$$

(29)

The simplexes created from the parameters with a length dimension are as follows:

$${\pi }_4=\frac{S_{\Sigma }}{S_{1,\Sigma ,\mathrm{c}}}$$

(30)

$${\pi }_5=\frac{S_{\mathrm{e}}}{S_{1,\Sigma ,\mathrm{c}}}$$

(31)

$${\pi }_6=\frac{S_{\mathrm{a},\Sigma }}{S_{1,\Sigma ,\mathrm{c}}}$$

(32)

$${\pi }_7=\frac{L}{d_1}$$

(33)

The fouling layer is affected by the flow surface area; therefore, the following simplex was created:

$${\pi }_8=\frac{d_1^2}{{{(d}_1-2·h_f)}^2}$$

(34)

The simplex created from the gas and air density parameters is as follows:

$${\pi }_9=\frac{{\rho }_{\mathrm{N}\mathrm{G}}}{{\rho }_{\mathrm{a}}}$$

(35)

The simplex for the thermal conductivity coefficients for the gas and for the fouling layer is as follows:

$${\pi }_{10}=\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{N}\mathrm{G}}}$$

(36)

With regard to the created simplexes, in the set of Equations (22)–(25), x₂ = x₃ = x₅ = x₈ = x₉ = x₁₀ = x₁₁ = x₁₂ = x₁₄ = x₁₆ = 0; therefore, the resulting equations are as follows:

$${ x}_1+x_{13}{+x}_{15}=0$$

(37)

$$2x_1+2x_6+x_7-{3x}_{13}+x_{15}=0$$

(38)

$$x_1+x_{15}=0$$

(39)

$$-x_1{+x}_4-x_{15}=0$$

(40)

The total number of searched criteria is l = 12 and the number of created simplexes is 10. The two missing complex criteria had to be identified based on two independent solutions of the set of Equations (37)–(40).

The set of Equations (37)–(40) contains six unknowns; therefore, they must be solved by determining two of those unknowns and calculating the remaining four.

In the first option, the determined unknowns were x₁ = 1 and x₆ = 0 and the calculated parameters were the following: x₄ = 0; x₇ = −1; x₁₃ = 0; and x₁₅ = −1. The identified criterion was as follows:

$${\pi }_{11}=\frac{K_{\mathrm{N}\mathrm{G}}}{{\lambda }_{\mathrm{N}\mathrm{G}}·d_1 }$$

(41)

In the second option, the determined unknowns were x₆ = 1 and x₁ = 0 and the calculated parameters were the following: x₄ = 0; x₇ = −2; x₁₃ = 0; and x₁₅ = 0. The identified criterion was as follows:

$${\pi }_{12}=\frac{S_{1,\Sigma ,\mathrm{c}}}{d_1^2 }$$

(42)

Since the criterial equation was to be derived for a particular cooler, the criteria π₄, π₅, π₆, π₇, and π₁₂ were constant. Therefore, they were pooled into a single new criterion π₀.

As a rule, similarity criteria may be transformed into other criteria through multiplication, division, exponentiation with a constant, or multiplication by a constant [21]. That rule was applied to obtain the following new criterion from the constant criteria:

$${\pi }_0=\frac{{\pi }_4·{\pi }_7}{ {\pi }_5·{\pi }_6·{\pi }_{12}}$$

(43)

The parameter that was being identified was the thickness of the fouling layer in the cooler tubes h_f, and since it fell within the π₈ criterion, it was expressed as a function of other criteria, with the following form:

$${\pi }_8=\psi ({\pi }_0,{\pi }_1,{\pi }_2,{\pi }_3,{\pi }_9,{\pi }_{10},{\pi }_{11})$$

(44)

The correlation between the dimensionless arguments in Equation (44) is exponential; therefore, the resulting equation is as follows:

$${\pi }_8=C· {\pi }_0^{z_0}·{\pi }_1^{z_1}{·\pi }_2^{z_2}{·\pi }_3^{z_3}{·\pi }_9^{z_9}{·\pi }_{10}^{z_{10}}{·\pi }_{11}^{z_{11}}$$

(45)

In a logarithmic scale, that correlation is linear and in the following form:

$${\mathrm{l}\mathrm{n}\pi }_8=\mathrm{l}\mathrm{n}C+{z_0·\ln {\pi }_0+z}_1·\ln {\pi }_1+z_2·\ln {\pi }_2+z_3·\ln {\pi }_3+z_9·\ln {\pi }_9+z_{10}·\ln {\pi }_{10}+z_{11}·\ln {\pi }_{11}$$

(46)

The individual criteria were calculated from the parameter values that were identified through iterations as described above. For the multiple regression analysis, data from 94 iterations were available. The C constant and the unknown exponents z_j were identified through multiple linear regression.

The coefficient of determination for the multiple linear regression was 0.9999. The regression sum of squares was 13.153, while the residual sum of squares was 0.005. The values of the C constant and the individual exponents are listed in

Table 3. C constant and the individual exponents.

Constant	Exponent
C	z0	z1	z2	z3	z9	z10	z11
1.3848 × 10−72	0.0000	0.017396	1.4880	−78.565	−83.971	173.23	0.0000

.

For the analysed cooler, the thickness of the fouling layer in the cooler tubes was subjected to the following criterial Equation (47), obtained by breaking down Equation (45):

$$\begin{array}{cc}\hfill \frac{d_1^2}{{{(d}_1-2{·h}_{\mathrm{f}})}^2}=& 1.3848· {10}^{-72}·{\left(\frac{K_{\mathrm{N}\mathrm{G}}}{K_{\mathrm{a}}}\right)}^{0.017396}· {\left(\frac{T_{1,\mathrm{N}\mathrm{G}}}{{\mathsf{\Delta }T}_{\mathrm{N}\mathrm{G}}}\right)}^{1.4880}·{\left(\frac{T_{2,\mathrm{N}\mathrm{G}}}{T_{1,\mathrm{a}}}\right)}^{-78.565}·\hfill \\ & ·{\left(\frac{{\rho }_{\mathrm{N}\mathrm{G}}}{{\rho }_{\mathrm{a}}}\right)}^{-83.971}·{\left(\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{N}\mathrm{G}}}\right)}^{173.23} \hfill \end{array}$$

(47)

Criteria with a zero exponent were not specified in Equation (47). In the linear regression, their impact was automatically included in the value of the C constant.

The parameter that was being identified—the fouling layer thickness h_f—is on the left side of criterial Equation (47) and is affected by the internal diameter of the cooler tube d₁. An adjustment of this criterion leads to solving a quadratic equation. Only one of the identified roots of the equation had a physical meaning—the one with a lower value. The values of the fouling layer thickness obtained from the model and from the analytical solution exhibited very good concordance (Figure 4).

The correlation may be described by a regression line with a slope approaching 1, more specifically 0.9987, with a reliability value (a squared correlation index) R² = 0.9999. The standard deviation of the difference was 0.0331, i.e., 3.2%.

The results of both solutions were also subjected to a pairwise t-test at a significance level α = 0.05. A precondition for using a pairwise t-test was meeting the assumption of normality of the difference scores, which was verified in a Shapiro–Wilk test of normality. The value of the test criterion t was identified using the following equation:

$$t=\frac{\left|\overline{{∆h}_{\mathrm{f}}}\right|·\sqrt{m-1}}{s_{∆}}$$

(48)

wherein m is the number of values (1) and s_Δ is the standard deviation (mm).

Parameter $\overline{{∆h}_{\mathrm{f}}}$ is the average value of the difference Δh_f,analyt,i − Δh_f,mod,i, calculated using the following equation:

$$\overline{{∆h}_{\mathrm{f}}}=\frac{\sum \left(h_{\mathrm{f},\mathrm{a}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{y}\mathrm{t},i}-h_{\mathrm{f},\mathrm{m}\mathrm{o}\mathrm{d},i}\right)}{m}$$

(49)

For the analysed data sets, m = 94 and $\overline{{∆h}_{\mathrm{f}}}=0.000102$ mm; therefore, the value of the test criterion t was 0.030. The critical value t_cr at a significance level of 0.05 was t_0.05(94−1) = 1.986. Since t < t_cr, it may be stated that both methodologies provided identical results.

2.3. Multiple Regression Analysis

In a general case, the derived Equation (47) describes the fouling layer thickness with sufficient accuracy. However, the model is not very convenient for common purposes as it includes 12 independent parameters. Therefore, the investigation was carried out with the aim of exploring how to reduce the number of relevant parameters. In any cooler, the individual criteria of Equation (47) contain certain parameters that are of a constant value. In the analysed cooler, those parameters included, for example, the thermal conductivity coefficient for the pipe material and the fouling layer, the internal and external diameters of the pipe, the pipe lengths, etc.

In this case, the fouling layer thickness was identified by applying multiple regression analysis, using the MinitabX 18 statistical software and the R package software 4.3.3 [22].

With the use of MinitabX, out of the group of relevant parameters, the following six parameters were selected as those with the most significant impact on the fouling layer thickness: Q_m,a; T_2,NG; T_1,a; ρ_NG; ρ_a; and c_p,NG. Multiple linear regression was applied and the following model equation was created:

$$\begin{gathered} h_{\mathrm{f}}=-5784.25-0.00491·Q_{m,\mathrm{a}}+1.32825·T_{2,\mathrm{N}\mathrm{G}}{+0.45745·T}_{1,\mathrm{a}}+6.90863·{\varrho }_{\mathrm{N}\mathrm{G}}+113.457·{\varrho }_{\mathrm{a}}+ \\ +1.70938·c_{p,\mathrm{N}\mathrm{G}} \end{gathered}$$

(50)

Gas density ρ_NG was used in Equation (50) for the mean temperature value (T_1,NG − T_2,NG)/2, and it was identified using Equation (51):

$${\rho }_{\mathrm{N}\mathrm{G}}=1.4563·{10}^{-3}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{N}\mathrm{G}}^2-3.6377·{10}^{-1}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}.\mathrm{N}\mathrm{G}}+63.00$$

(51)

Air density ρ_a was calculated using Equation (52):

$${\rho }_{\mathrm{a}}=1.313·{10}^{-5}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{a}}^2+4.628·{10}^{-3}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}.\mathrm{a}}-1.276$$

(52)

wherein $t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}.\mathrm{a}}=t_{1,\mathrm{a}}+{\mathsf{\Delta }T}_{\mathrm{N}\mathrm{G}}$.

The specific thermal capacity of gas c_p,NG was calculated using Equation (53):

$$\begin{gathered} c_{p,\mathrm{N}\mathrm{G}}=1.1696·{10}^{-9}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{N}\mathrm{G}}^4-6.3457·{10}^{-7}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{N}\mathrm{G}}^3+1.3139·{10}^{-4}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{N}\mathrm{G}}^2- \\ -1.0701 ·{10}^{-2}·t_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n},\mathrm{N}\mathrm{G}}+2.9979 \end{gathered}$$

(53)

Equation (50) represents the model of a correlation between the fouling layer thickness and the selected input parameters. The regression model parameters were identified by applying the method of least squares, which is sensitive to outliers. Therefore, the input data was also assessed in terms of outliers and influential values, which were excluded from the model.

The statistical significance of the regression model, or of the regression model parameters, was verified in the tests of statistical significance at a significance level α = 0.05. As a rule, if the p-value is lower than the significance level α, then the null hypothesis is rejected in favour of the alternative hypothesis. If the p-value equals or is higher than the selected significance level α, then the null hypothesis is not rejected [23]. The results of the testing showed that the regression model (50) was statistically significant (p-value = 3 × 10⁻¹²² < α) and that all of the analysed parameters were also statistically significant (p-value < α). The coefficient of determination for the multiple linear regression was 0.9985; this means that as much as 99.85% of the fouling layer thickness variability may be explained by the proposed regression model with the given parameters.

Heteroscedasticity of the regression model was tested with the use of the Breusch–Pagan test, in which the null hypothesis assumes homoscedasticity. The results (p-value = 0.199 > α) indicated that the regression model did not exhibit heteroscedasticity.

The correlation between the fouling layer thickness obtained from the model created by applying the multiple regression and that obtained analytically may be described by a regression line (Figure 5) with a slope approaching 1, in particular 0.9996, at a reliability value R² = 0.9985.

The values of the fouling layer thickness obtained from the model and from the analytical solution were compared in a pairwise t-test at a significance level α = 0.05. The Shapiro–Wilk test of normality was used to verify data normality. The results of the pairwise t-test (p-value = 0.999 > α) indicated that the two methodologies provided comparable results. The tests were performed using the R package software [24].

3. Conclusions

Based on the performed analyses, it may be stated that the application of the balance method and dimensional analysis to the processes that are affected by a relatively large number of relevant parameters is not very beneficial. It is more appropriate to describe such processes with the use of a regression equation constructed from selected parameters, and not from similarity criteria.

The constructed regression Equation (50) is simple and easy to use for a double crossflow cooler of natural gas. It is therefore unnecessary to put the cooler out of service and dismantle the tubes inside the cooler in order to examine the extent of the fouling layer accumulated on the surface inside the tubes. Such interruptions in the cooler operation cause severe outages in the cooling process and result in a considerable financial loss associated with an interruption in the transportation of gas, as well as the dismantling and re-assembly of the cooler tubes, because there are thousands of tubes in a cooler.

Equation (50) was validated during the process of cleaning the analysed cooler. Prior to the cleaning, the gas temperature at the outlet from the cooler was 15 K lower than the temperature at the inlet. The temperature of the surrounding air was 20 °C. Based on the volume of the fouling matter that was pushed out of all the pipes in the cooler, the average thickness of the fouling layer on the inner pipe surface was calculated as 2.22 mm. The thickness of the fouling layer calculated using Equation (50) was 2.15 mm. The difference between the two values represents 3.1%.

Equation (50) may be simply programmed, and in case of a change in the gas output temperature and the ambient temperature, it facilitates identifying the current thickness of the fouling layer on the internal surface of the cooler tubes without the need to put the entire cooler out of service.

The presented methodology may be applied at any compressor station in Slovakia, as well as any other country, where natural gas is cooled by coolers of the same design as the design of the cooler subjected to the analysis described in this article.