Rapid Prediction of the In Situ Pyrolysis Performance of Tar-Rich Coal Using the POD Method

Abstract

In this paper, a POD reduced-order interpolation model for solving the in situ pyrolysis process of tar-rich coal is employed to predict the flow and heat transfer performance in the porous media region so as to save computational resources and realize fast calculations. Numerical simulation using the finite volume method (FVM) is firstly used to obtain sample data, based on the samples through the primary function and spectral coefficients of the solutions. The physical field information and parameter distribution under different conditions of inlet temperature, inlet velocity and permeability are predicted. The results are compared with those of FVM to verify the accuracy of the calculated results. The relative mean deviation (RME) of the results of the POD prediction of each parameter for each working condition was synthesized to be no more than 5%. The performance of in situ pyrolysis of tar-rich coal is then investigated, and the oil and gas production are predicted. As the inlet velocity increases from 0.3 m/s to 0.9 m/s, the fraction of high-quality oil and gas production reaches 0.47 and then decreases to 0.38. Increasing the inlet temperature and permeability has a negative effect on the fraction of high-quality hydrocarbon production, after which the quality fraction of high-quality oil and gas dropped sharply to about 0.22. Porosity has a positive impact on the oil and gas production. When the porosity reaches 0.3, the quality fraction of high-quality oil and gas can reach 0.27.

1. Introduction

Tar-rich coal refers to coal with a dry tar base yield of 7–12%. Its biggest feature is that coal is rich in hydrogen-rich structure, which can generate oil and gas by pyrolysis, and has a high tar yield. It is a special coal resource that integrates coal and oil properties [1]. In situ pyrolysis of tar-rich coal is a new production method of directly heating the original coal reservoir through heat carriers to promote the underground pyrolysis of tar-rich coal [2]. During the chemical process of pyrolysis, the heat is transferred from the heat carrier injected through the inlet well to the coal reservoir, and the reaction occurs after the local temperature of the coal is higher than the pyrolysis temperature. Since the coal reservoir is porous and there exist massive pores and cracks inside the coal seam, both the heat carrier such as nitrogen or water vapor and the pyrolysis products such as the gaseous tar, methane, carbon monoxide, carbon dioxide and hydrogen would diffuse through the pores and cracks to the production well [3]. Some research has focused on coal seam gas [4,5,6,7]. For example, Wang et al. [4,5] studied the influence of negative pressure on the leakage during the extraction of coalbed methane and found that the active support method using delayed expansion material can effectively reduce the formation of crack micropores and reduce stress concentration. The technology of sealing holes with three plugging and two filling solid and liquid materials is put forward to realize continuous high-concentration gas extraction [6]. Zhu et al. [7] developed a simulation model to evaluate the effect of penetrating fractures on lateral water flow and to study the effect of permeability on gas storage and production in coal seams. The above literature illustrates the risk of leakage in coal seam gas extraction. At the same time, the temperature of the coal reservoir varies with the heating by the heat carriers and the endothermic reaction. Thus, in situ pyrolysis is a rather complicated multi-physics process including heat and mass transfer as well as chemical reactions inside the porous media.

Due to the high cost, high risk, harsh conditions and difficult underground measurement of in situ pyrolysis of tar-rich coal, it is difficult to meet the demands of practical applications. Therefore, it is significant to predict the performance of in situ pyrolysis of tar-rich coal under unknown working conditions for engineering applications. Numerical simulation based on the finite volume method (FVM) can obtain the heat and mass transfer characteristic, but it is time-consuming and consumes large computational resources. Proper orthogonal decomposition (POD) is a fast, efficient and reliable method to predict unknown working conditions based on known data [8]. It obtains the key features from a series of samples, which ensures the computational accuracy and saves the computational resources at the same time. The POD algorithm includes the interpolation and the projection method. The difference between the two methods is reflected in the different solving paths of spectral coefficients. POD interpolation relies on the number of variables, while POD projection relies on the governing equation for solving [9]. For example, the dominant coherent structure of cavitation flow is studied using the POD and DMD methods by Liu et al. [10]. Lu et al. [11] applied the POD method to double rotor bearing downscaling and investigated the frequency characteristics of the deviation response of a double rotor bearing coupling. Sun et al. [12] predicted the bending shock dominant flow field at the inlet of a punching machine under different attack angles and free flow Mach number using the POD method. In addition, the POD method has also been applied to problems such as fin-tip vortex [13], VCT (vibration correlation technique) [14], rough plate transient flow [15] and other problems. Then, the POD-Galerkin method was proposed and applied in some research to improve the calculation accuracy [16,17,18,19]. In addition, combining with other interpolation methods or models, the POD method can improve the model and obtain more accurate results [20,21,22,23,24].

Compared with the POD projection method, the POD interpolation method does not require the derivation of differential equations. Therefore, the POD interpolation method has higher applicability and efficiency than the POD projection method. Moreover, the POD interpolation method has a simpler operation process, and the deviation is acceptable in most cases. For in situ pyrolysis of tar-rich coal, numerical simulation of the flow heat transfer process in porous media region is required. The POD method has also been reported to solve the flow heat transfer problem in the porous media region. For example, Li et al. [25,26] proposed a new POD-ROM to deal with the flow and heat transfer in the fractured region. The global order reduction model of steady-state flow of the fractured region is established by using orthogonal decomposition and the Galerkin projection method. Selimefendigil et al. [27] studied the convection drying performance of porous wet objects under channel flow and multi-impact jet structure and reduced the calculation cost greatly using the POD method.

In this paper, a POD reduced-order interpolation model for solving the in situ pyrolysis process of tar-rich coal is employed to predict the flow and heat transfer process in the porous media region so as to save computational resources and realize fast calculations. Numerical simulation using the finite volume method (FVM) is firstly used to obtain sample data, based on the samples through the primary function and spectral coefficients of the solutions. Then the POD prediction data are calculated, and the results are compared and analyzed to determine the prediction accuracy. The performance of in situ pyrolysis of tar-rich coal is then investigated, and the oil and gas production are predicted.

2. POD Reduced-Order Model

The implementation process of POD reduced-order model [8] is shown in Figure 1. The first sampling step is to build a sample matrix based on the data obtained from FVM. Then, the primary function is calculated based on the Snapshot method [28]. The spatial discrete points are much larger than the temporal discrete points. The spectral coefficients are then calculated using interpolation. Finally, the approximate solution space is obtained using linear reconstruction of the principal functions and spectral coefficients.

The sampling data is represented by the matrix S. As shown in Equation (1), each column stores only one sample corresponding to one moment of a particular set of parameters.

$$\mathrm{S}=\left[\begin{array}{cccccccccc}{p}_{1,1}^{1\left(t_0\right)}& p_{1,1}^{1\left(t_1\right)}& \cdots & p_{1,1}^{1\left(t_F\right)}& p_{1,1}^{2\left(t_0\right)}& \cdots & p_{1,j}^{2\left(t_F\right)}& p_{1,1}^{13\left(t_0\right)}& \cdots & p_{1,j}^{13\left(t_F\right)}\\ ⋮& ⋮& & ⋮& ⋮& & ⋮& ⋮& & ⋮\\ p_{i,1}^{1\left(t_0\right)}& p_{i,1}^{1\left(t_1\right)}& \cdots & p_{i,1}^{1\left(t_F\right)}& p_{i,1}^{2\left(t_0\right)}& \cdots & p_{i,1}^{2\left(t_F\right)}& p_{i,1}^{13\left(t_0\right)}& \cdots & p_{i,j}^{13\left(t_F\right)}\\ p_{1,2}^{1\left(t_0\right)}& p_{1,2}^{1\left(t_1\right)}& \cdots & p_{1,2}^{1\left(t_F\right)}& p_{1,2}^{2\left(t_0\right)}& \cdots & p_{1,2}^{2\left(t_F\right)}& p_{1,2}^{13\left(t_0\right)}& \cdots & p_{1,2}^{13\left(t_F\right)}\\ ⋮& ⋮& & ⋮& ⋮& & ⋮& ⋮& & ⋮\\ p_{i,2}^{1\left(t_0\right)}& p_{i,2}^{1\left(t_1\right)}& \cdots & p_{i,2}^{1\left(t_F\right)}& p_{i,2}^{2\left(t_0\right)}& \cdots & p_{i,2}^{2\left(t_F\right)}& p_{i,2}^{13\left(t_0\right)}& \cdots & p_{i,2}^{13\left(t_0\right)}\\ ⋮& ⋮& & ⋮& ⋮& & ⋮& ⋮& & ⋮\\ p_{i,j}^{1\left(t_0\right)}& p_{i,j}^{1\left(t_1\right)}& \cdots & p_{i,j}^{1\left(t_F\right)}& p_{i,j}^{2\left(t_0\right)}& \cdots & p_{i,j}^{2\left(t_F\right)}& p_{i,j}^{13\left(t_0\right)}& \cdots & p_{i,j}^{13\left(t_0\right)}\end{array}\right]$$

(1)

where $p_{i,j}^{number\left(t_k\right)}$ denotes the parameter that specifies a point in space and time, with the subscript denoting the spatial point at a specific operating condition. The superscript number indicates the number under particular operating conditions, and the superscript $\left(t_k\right)$ denotes a specific point in time. From left to right corresponds to the initial time of the sample to the final time, followed by arranging the sample data of other conditions. All the sample data of conditions representing different processes are given in matrix S. Based on the matrix S, some unknown conditions could also be predicted.

The working conditions and parameters for the two models obtained using FVM numerical simulations are shown in

Table 1. Working conditions of fractured zone of in situ pyrolysis.

Number	Inlet Flow Velocity	Temperature	Injection Mode
1	0.1 m/s	973 K	Gas injection between cracks
2	0.3 m/s	973 K	Gas injection between cracks
3	0.5 m/s	973 K	Gas injection between cracks
4	0.7 m/s	973 K	Gas injection between cracks
5	0.9 m/s	973 K	Gas injection between cracks
6	0.5 m/s	873 K	Gas injection between cracks
7	0.5 m/s	923 K	Gas injection between cracks
8	0.5 m/s	973 K	Gas injection between cracks
9	0.5 m/s	1023 K	Gas injection between cracks
10	0.5 m/s	1073 K	Gas injection between cracks

and

Table 2. Working conditions of homogeneous porous region of in situ pyrolysis.

Number	Permeability	Porosity	Inlet Pressure	Temperature
1	0.1 mD	0.26	200 kPa	973 K
2	0.2 mD	0.26	200 kPa	973 K
3	0.5 mD	0.26	200 kPa	973 K
4	2 mD	0.26	200 kPa	973 K
5	5 mD	0.26	200 kPa	973 K

. The FVM simulation is conducted according to the previous paper [3] and is not introduced in detail here. According to the working conditions listed in the tables, the sample data are arranged sequentially in order to obtain the sample matrix. Each sample data represents the simulation data under the working condition corresponding to a specific parameter.

The Snapshot method [29] is used to construct the following inner product matrix C to obtain the primary function.

$$C=S{S}^T$$

(2)

The eigenvalues λ and corresponding eigenvectors V are obtained by singular value decomposition as shown in Equation (3). Thus, the primary function can be calculated from Equation (4).

$$CV=\lambda V$$

(3)

$$\phi =\frac{SV}{\Vert SV\Vert }$$

(4)

Subsequently, the number of primary functions used in POD is determined by calculating the energy contribution and cumulative energy contribution at different numbers of primary functions to achieve higher prediction accuracy while improving computational efficiency.

The energy contribution rate and cumulative energy contribution rate are shown in Equations (5) and (6), respectively.

$${\xi }_i=\frac{{\lambda }_i}{{\sum }_{k=1}^N{\lambda }_k}$$

(5)

$$\epsilon =\sum _{i=1}^M{\xi }_i=\frac{{\sum }_{k=1}^M{\lambda }_k}{{\sum }_{k=1}^N{\lambda }_k},M\le N$$

(6)

where N and M are the number of total primary functions and selected primary functions, respectively; ${\lambda }_k$ is sorted in descending order; ${\xi }_i$ is the energy contribution of the $i$th primary function; $\epsilon$ represents the cumulative energy contribution of the POD primary functions. In practice, when $\epsilon$ is closer to 1, the results are more accurate.

Using the POD interpolation method, the spectral coefficients can be realized. Assume $f\left(x,y,t,p_n\right)$ $\left(n =1, \dots , N\right)$ represents the physical field of parameter $p$ and $n$ denotes the number of parameters $p$. A set of primary functions ${\phi }_n(x,y)$ can be obtained by performing POD prediction on these N samples. The corresponding spectral coefficients under the design parameters can be obtained from Equation (7).

$$c\left(p_n\right)=(f\left(x,y,t,p_n\right) {\phi }_n(x,y\left)\right)$$

(7)

The spectral coefficients of unknown circumstance can be calculated by Equation (8).

$$c(p)=c(p_n)+[c(p_{n+1})-c(p_n)]\frac{p -p_n}{{p_{n+1} -p}_n}$$

(8)

$$f\left(x,y,t,p_n\right)=\sum _{k=1}^{M}c(p_n){\phi }_n(x,y)$$

(9)

Finally, any parameter of the demand solution and its distribution can be obtained through Equation (9).

In this study, N takes the value of 5. The five sampling points include the inlet temperature, outlet temperature, and three other sampling points. The locations of the remaining three sampling points are shown in Figure 2, representing the temperature near the center and edge at the corresponding time of three months, respectively.

After the POD calculation, the deviation relative to the FVM calculation needs to be analyzed. The absolute deviation is calculated by comparing the data from POD prediction and FVM numerical simulation and under the same operating conditions.

$$E=\frac{{\sum }_{i=1}^n(p_{POD}-{ p}_{FVM})}{n}$$

(10)

where $p_{POD}$ and $p_{FVM}$ are the parameter values for POD prediction and FVM simulation, respectively.

Relative Mean Error (RME) can be calculated by Equation (11).

$$RME=\left|\frac{E}{p}\right|\times 100\%$$

(11)

In order to reveal the relation between POD prediction parameters and FVM simulation parameters, Equation (12) was used with reference to Fang et al. [29].

$$R=\frac{Cov(p_{POD},p_{FVM})}{\sqrt{Var\left(p_{POD}\right)Var\left(p_{FVM}\right)}}$$

(12)

where $Cov(p_{POD},p_{FVM})$ is the covariance of $p_{POD}$ and $p_{FVM}$, $Var\left(p_{POD}\right)$ is the variance of $p_{POD}$, and $Var\left(p_{FVM}\right)$ is the variance of $p_{FVM}$.

3. Results and Discussion

In this section, the POD reduced-order model is defined to solve the in situ pyrolysis process of tar-rich coal and predict the physical field information and parameter distributions under different working conditions of inlet temperature, inlet flow velocity and permeability. The results obtained by the POD prediction are compared with the FVM results, and the accuracy is verified by calculating the correlation coefficient and the average relative deviation pair to ensure the reliability and accuracy of the POD calculation. Then, the oil and gas production is predicted based on the POD calculation.

3.1. Comparative Analysis of Inlet Velocity Prediction Results

The POD prediction is performed at the inlet flow velocity of 0.5 m/s and the inlet temperature of 973 K (corresponding to Case 3 in Table 1). The initial temperature and pressure of the coal seam are 293 K and 0 Pa, respectively. In the following sections, the same initial state is kept. Firstly, the effect of primary function is investigated, and the results are shown in Figure 3. The cumulative energy contribution gradually improves from 10.70% to 100% as the number of selected primary functions increases. The energy contribution reaches 98.77% using four primary functions for which the calculation accuracy and calculation efficiency are both acceptable. Thus, four primary functions are adopted in the following cases.

The temperature distribution comparison between the POD and the FVM results are shown in Figure 4. From the figure, we can see that the temperature distribution obtained from the POD prediction is approximately the same as that from the FVM. At the junction of the high-temperature and low-temperature zones, where the temperature changes are the most drastic, there is a small difference between the calculation results of POD and FVM, which is due to the intense heat transfer, the fast rate of pyrolysis reaction and the large temperature change over time in this region, and the interpolation difference generated by the POD as well as the interpolation errors in the process of the subsequent cloud mapping leads to the calculation difference in this place. The distribution of the POD results is relatively uniform and can save computational resources while reflecting the trend of temperature change with small deviation, which is of practical significance.

The temperature profiles at y = 5 m and x = 5 m calculated by POD and FVM are further compared at 3 months, as shown in Figure 5. The results show that at y = 5 m, the temperature trends calculated by POD and FVM are basically the same, but due to the temperature change in the middle section of y = 5 m and the high chemical reaction rate, the difference in the results obtained by POD prediction is relatively large, and the average relative deviation of the overall temperature is 3.46%, which is a small deviation; while at x = 5 m, it is known from its temperature field that the central region is close to the inlet temperature 973 K, while the nitrogen passing through the edge is relatively less, the heat transfer and reaction are weaker, and the temperature is lower. Thus, the difference between the POD calculation results and the FVM calculation results becomes larger, and the average relative deviation of the overall temperature is 0.53%, which is smaller than that for y = 5 m. The average relative deviation is 3.46%.

The cumulative energy contribution at different number of primary functions for the oil and gas mass fraction is shown in Figure 6. Four primary functions are used, with the cumulative energy contribution being 97.88%.

The distribution of oil and gas mass fraction obtained by FVM and POD at the inlet velocity of 0.5 m/s is shown in Figure 7. Due to the fast diffusion of gas and the influence of other factors such as reaction and other perturbations, there are some differences between the POD prediction results and the FVM calculation results. At 3 months, the region near the inlet has been fully reacted, the main reaction area is distributed on the two sides above and below the outlet, and the reaction is relatively mild, thus, the difference between the POD results and the FVM results is smaller. The overall mass fraction distribution is similar, and the accuracy of POD prediction is acceptable in these cases.

3.2. Comparative Analysis of Inlet Temperature Prediction Results

The results comparison at the inlet temperature of 973 K and the inlet flow velocity of 0.5 m/s (corresponding to Case 8 in Table 1) are shown in Figure 8. It can be seen that the temperature field distributions of POD and FVM turn out to be approximately the same. The POD calculation can reflect the temperature change trend, and the deviation is small when compared with the FVM calculation. At the junction of the high-temperature zone and low-temperature zone, the temperature change is the most violent, which is due to the violent heat transfer in this region, the fast rate of pyrolysis reaction and the large temperature change over time. There is a difference between the calculation results of POD and FVM, which is caused by the interpolation errors generated by POD and cloud mapping.

Similarly, the distribution results of oil and gas mass fraction can be calculated by building sample matrices based on the data at different inlet temperatures. The differences between the POD prediction results and the FVM results are also small. The distribution results have the same trend over time calculated by the two methods. When three POD primary functions are selected, the cumulative energy contribution reaches 98.50% as shown in Figure 9, and the computational accuracy and computational efficiency are high.

As shown in Figure 10, the results for oil and gas mass fraction at the inlet temperature of 973 K by two methods are basically the same, and the deviation of the results is very small with a maximum deviation of 0.01%. At this time, only two primary functions need to be selected, corresponding to the cumulative energy contribution having reached 94%, with high computational accuracy and computational efficiency, as shown in Figure 11.

3.3. Comparative Analysis of Permeability Prediction Results

When the inlet pressure is 200 kPa, the inlet temperature is 973 K, the permeability is 0.5 mD and the porosity is 0.26 (corresponding to Case 3 in Table 2), POD prediction calculations are performed and compared with the FVM calculations.

From Figure 12, we can see that the temperature distribution results of the two methods are similar, and the POD prediction accuracy is higher. At the junction of the high- and low-temperature zones where the temperature change is drastic, the deviation of the calculation is also higher.

The comparison of POD and FVM results at the permeability of 0.5 mD is shown in Figure 13. The results show that the trends are basically the same with a small deviation.

From the above cases, we can obtain that the average computation time of in situ pyrolysis of tar-rich coal for POD is about 5.25 s, while the average computation time for FVM is estimated to be 54 h. In terms of computation time, POD computation leads to an increase in computational efficiency by a factor of about 3.88 × 10⁴. The field data of unknown conditions can be obtained quickly by the POD method. According to the existing data of simulation, the data of the whole field can be given by the POD algorithm quickly with high accuracy, such as the exit steam mass fraction, average temperature, etc.

3.4. Deviation Analysis

The correlation coefficients and mean relative deviations of the temperature calculated by FVM and POD at an inlet velocity of 0.5 m/s, temperature of 973 K and permeability of 0.5 mD are shown in Figure 14a. The results show that the POD predicted temperature has a strong correlation with the FVM simulation under this working condition, and the correlation coefficient is greater than 0.75 in 1–5 months. As shown in Figure 14b, the average relative deviation stays around 4.73% during 1–5 months, the correlation coefficient is greater than 0.88 within 1–5 months under this condition and the prediction of temperature by POD is strongly correlated with the FVM calculation. The maximum average relative deviation is less than 2.2% in the first 5 months, indicating better accuracy and reliability of the POD prediction. In Figure 14c, the correlation coefficient is greater than 0.97 during the period of 1–5 months; the relative mean deviation is maintained at about 4.70% during the period of 1–5 months; the relative mean deviation (RME) of the POD prediction for each parameter of each working condition is no more than 5%.

3.5. POD Prediction of Oil and Gas Quality in the Production

Due to the high cost, high risk, harsh conditions and difficult underground measurement of in situ pyrolysis of tar-rich coal, it is difficult to meet the demands of practical applications. Proper orthogonal decomposition (POD) is a fast, efficient and reliable method to predict unknown working conditions. Based on this method, the heat and mass transfer process of pyrolysis can be efficiently predicted. In addition, we also pay attention to the obtained tar during pyrolysis. The POD method can also predict the quality of oil products to obtain high-quality oil products.

In this section, POD is employed to predict the performance of pyrolysis for more unsimulated cases. In order to analyze the effect on in situ pyrolysis of tar-rich coal comprehensively, the fraction of high-quality oil and gas production is defined to reflect the quality of oil and gas produced by in situ pyrolysis, as given in Equation (13). The fraction of high-quality oil and gas production $z$ represents the ratio of total amount of oil and gas production in the temperature range of 673–873 K to the whole production during the process. The higher the fraction of high-quality oil and gas production in this temperature stage, the larger the economic gain obtained.

$$z=\frac{{\int }_{673}^{873}\frac{d{O}_{tar}}{dT}dT}{O_{tar}}$$

(13)

As shown in Figure 15a, with the increase in inlet flow velocity, the fraction of high-quality oil and gas production fraction reaches 0.47 and then decreases to 0.38. The increasing inlet flow velocity strengthens the heat transfer at the beginning. Subsequently, with further increase in inlet flow velocity, the temperature in some areas rises rapidly and the uniformity of heat transfer becomes worse. Thus, the yield fraction of high-quality oil and gas is reduced. The maximum fraction of high-quality oil and gas production is obtained under the condition that the inlet flow velocity is 0.5 m/s. Figure 15b shows the variation of oil and gas production fraction with different inlet temperature. As the inlet temperature rises from 873 K to 1073 K, the high-quality oil and gas production fraction gradually decreases. This is because the heat exchange between nitrogen and coal seam gets more rapid and intense at higher inlet temperature and the oil and gas output in the main reaction stage decreases. As shown in Figure 15c, when the porosity reaches 0.30, the high-quality oil and gas production fraction gradually increases to 0.27, and the diffusion rate of nitrogen increases and the flow and heat transfer are strengthened, which is conducive to producing more hydrocarbons in the same temperature variation range. Therefore, the high-quality hydrocarbon production fraction can be obtained. As for the permeability, it has a negative effect on the improvement of high-quality oil and gas production fraction. When the permeability changes from 0.1 to 5 mD, the fraction of high-quality oil and gas production decreases from 0.46 to 0.04, as shown in Figure 15d. The convective heat transfer is enhanced, and the temperature of the coal reservoir rises rapidly to the inlet temperature. The percentage of oil production in the temperature range from 673 K to 873 K also decreases.

4. Conclusions

On the basis of multi-physics coupled FVM numerical simulation of in situ pyrolysis of tar-rich coal, the sample matrix is constructed and a POD reduced-order model is established in this paper to solve the in situ pyrolysis process of tar-rich coal. The physical field information and parameter distributions corresponding to inlet temperature, inlet velocity and permeability are predicted. The main conclusions are as follows:

• The POD calculation can improve the calculation efficiency by about 3.88 × 10⁴ times with a relative mean deviation (RME) of no more than 5%.
• When the inlet flow velocity is 0.5 m/s, the high-quality oil and gas production fraction is the largest, and the oil and gas production efficiency is the highest.
• The fraction of high-quality hydrocarbon production decreases from 0.62 to 0.21 when the temperature rises from 873 K to 1073 K and changes from 0.14 to 0.27 with the porosity increasing from 0.18 to 0.30.
• The increase in permeability significantly enhances the convective heat transfer and the temperature of the coal reservoir rises rapidly, which has a negative effect on the fraction of high-quality oil and gas production.