The Effective Thermal Conductivity of Unsaturated Porous Media Deduced by Pore-Scale SPH Simulation

1 Introduction

Heat transfer in porous media is a common phenomenon during agricultural and industrial production (Moradi et al., 2016; Varnosfaderani et al., 2016; Cui et al., 2020). The effective thermal conductivity is a key problem for the heat-conducting properties of porous media (Jahangir et al., 2018a; Bai et al., 2021a; Fu et al., 2020), which involves the coupling effects of multi-physical fields. Changes in the thermal field result in changes in the moisture and stress fields of porous media (Das and Cleary, 2016; Jahangir et al., 2018b; Meng et al., 2020; Nishad and AI-Raoush, 2021; Yuan et al., 2021a). Therefore, it is important to determine the effective thermal conductivities of porous media and clarify how they relate to the influencing factors.

As a mixture, a porous medium is composed of solid components that form a solid framework and pore space that is filled by other components (Cheng et al., 2020; Hu et al., 2021; Yuan et al., 2021b). The random and disordered natures of these structures make it difficult to accurately quantify the effective thermal conductivities of porous media in theory. The effective thermal conductivity is a general parameter for characterizing porous media and is mainly affected by the components of the medium, more specifically, the distributions, proportions and properties of the components. When the solid, liquid and gas phases have a layered distribution and are connected in series, the thermal conductivity reaches its minimum; in contrast, when the layers are connected in parallel, the thermal conductivity reaches its maximum. When the layers are connected in series (superscript L) or in parallel (superscript U), the theoretical values of the effective thermal conductivity for a multiphase medium are

where λ is the thermal conductivity of the porous medium, the subscript n refers to different components, and $\varphi$ denotes the volume proportion of each component.

A disordered multiphase mixture can ideally be modeled as some components in series and some components in parallel. Based on this theory, the thermal conductivity values for multiphase media are within the abovementioned limits. For unsaturated porous media, the upper and lower limits of the effective thermal conductivity can be further refined. Since the thermal conductivity of water is much greater than that of air, the thermal conductivity of a porous medium increases with increasing water content. When in a saturated state, the thermal conductivity value of porous media reaches its upper limit; while for a completely dry state, the thermal conductivity reaches the lower limit.

Johansen (1975) proposed the concept of the normalized thermal conductivity λ_r, which is defined as

where λ_dry and λ_sat indicate dry and saturated conditions, respectively.

The value of λ_r varies between 0 and 1. By fitting experimental data from the literature, Johansen (1975) proposed a logarithmic relationship between λ_r and the degree of saturation S_r for several types of soil. By analyzing a large amount of literature data, Cote and Konrad (2005) developed the following empirical equation:

This equation can be transformed into

where α is an influence coefficient.

Through laboratory tests on quartz sand, Chen (2008) found that the logarithm of thermal conductivity for sand has a good linear relationship with porosity, and the slope is affected by the degree of saturation. This empirical model is

where $\phi$ is the porosity, b and c are both empirical parameters.

At present, studies on the thermal conductivity of soil have been mainly conducted at the macroscopic scale (Sass et al., 1971; Vincent and Paul., 2005; Bai et al., 2018). To achieve a deeper understanding of the laws governing the effective thermal conductivity, it is necessary to advance research to the pore scale (Cui et al., 2019; Zhou et al., 2022). This study employs a smoothed particle hydrodynamics (SPH) method due to its flexibility in pore-scale simulations and good computational accuracy in heat conduction calculations.

SPH is a kind of meshless particle algorithm of Lagrangian form. Initially proposed in 1977 to solve a problem in astrophysics, SPH has since been used in heat conduction. Cleary and Monaghan (1999) proposed an improved form of the standard SPH formula that ensures a heat flow balance on material interfaces, and the results showed that the formula is highly accurate. Jeong et al. (2003) improved the SPH method for solving second-order differential equations by decomposing the second-order partial differential equation into two first-order partial differential equations and then checked the accuracy by performing 1D and 2D heat conduction simulations. Jiang and Sousa (2007) employed the SPH method to predict the effective thermal conductivity in heterogeneous materials containing two or three different components. By combining the advantages of SPH and the finite difference method (FDM), Bai et al. (2018) proposed a well-established SPH-FDM method and analyzed the thermal processes of homogeneous media containing discontinuous interfaces and showed that this method can be successfully employed in parabolic problems.

The present work utilizes a pore-scale simulation approach to predict the thermal conductivities of unsaturated porous media using the SPH method and examines the relationship between the normalized thermal conductivity and saturation by using an empirical model. The manuscript is structured as follows: Section 2 introduces the reconstruction process of an unsaturated granular soil at the pore scale and the procedure for the simulation experiment; Section 3 verifies the accuracy of the simulation experiment through ideal media; and Section 4 presents simulation experiments in 3D granular media and discusses the variation in the normalized thermal conductivity with saturation. The conclusions drawn from the work are presented in Section 5.

2 Pore-Scale Simulation on Effective Thermal Conductivity

2.1 Smoothed Particle Hydrodynamics Method

The SPH method is a kind of meshless particle algorithm of Lagrangian form. The computational domain is discretized by separate sets of particles. The discretization of the field function consists of two main steps. First, the field function f at point x can be expressed in integral form:

where Ω, W, and h denote the integral domain, the kernel function and the smoothing length of the kernel function, respectively.

A frequently used kernel function, the cubic B-spline function, is adopted herein (Huang and Liu, 2020):

where R is the relative distance and is equal to r/h, where r is the particle distance. In 1D, 2D, and 3D spaces, the normalized values of the constant α_d are 1/h, 15/7πh², and 3/2πh³, respectively.

Next, the integral form of the field function is discretized into a series summation. The field function at particle i is approximated as

where x is the distance vector, m and ρ are the mass and density of particle j, respectively, the subscript j denotes a particle within the influence area of particle i, and N is the total number of particles.

However, the traditional SPH algorithm has several drawbacks when applied to complex engineering problems. To address these defects, many researchers have proposed improved calculation methods, which is a driving force for the further development of the SPH method (Huang et al., 2019; Rao and Bai, 2020). At the same time, this method is also being used in an increasingly broad range of applications, such as heat conduction problems at the pore scale and thermal consolidation problems in geotechnical materials (Bai et al., 2021b; Xue et al., 2021).

The differential equation for unsteady heat conduction within a three-dimensional isotropic medium without an internal heat source is (Carslaw and Jaeger, 1959).

where c is the specific heat.

The SPH approximation of Eq. 9 at point i can be written as (Sass et al., 1971; Bai et al., 2018).

where r_ij =|x_ij |.

The size of each time interval is suggested to be (Bai et al., 2018).

where κ is the thermal diffusion coefficient equal to λ/ρc.

The explicit Euler method is employed for the integration step, $\tilde{f}\left(t+\mathrm{\Delta }t\right)=f\left(t\right)+\mathrm{\Delta }t{\frac{df}{dt}|}_t$. SPH has different discretization schemes for different derivation processes (Price, 2012). Among these schemes, Eq. 10 is the most suitable for modeling thermal conduction. Its accuracy in 1D, 2D, and 3D spaces has been verified by extensive numerical tests (Cleary and Monaghan, 1999).

2.2 Generation of 3D Porous Media

The particle flow code (PFC) method is employed to produce a 3D granular medium. Cundall and Strack (1979) developed the PFC method from the discrete element method by defining the discrete element as a rigid sphere. Based on PFC3D 3.10 software (developed by ITASCA Inc., United States), there are five steps involved in generating a 3D dry granular medium: first, the software sets rigid walls that represent the boundaries of the sample; second, the software sets the number and size range of the spheres and randomly generates them in the calculation domain; then, gravity is applied to the spheres, which are allowed to fall freely; fourth, when the spheres have naturally accumulated, the software applies a plate with slow downward movement over the accumulated body to compact it to the predefined porosity; finally, the software outputs the radii of all the spheres and the coordinates of the centers of mass and then discretizes the medium into solid and gas particles (Figure 1A).

FIGURE 1. Solid and gas particles: (A) Restructured 3D porous medium (named SAP1); (B) discretization of SAP1 in the dry state (the green particles are gas particles, while the red particles are solid particles).

In Figure 1A, the length L of the cubic medium (named SAP1) is 0.008 m. It contains 70 spherical granules with identical radii of a = 1 mm. After spatial discretization, the medium (Figure 1B) contains a total of 64,000 particles with a particle spacing of 0.0002 m. The particles are divided into gas and solid particles. The porosity is set at 0.43. The calculation time is about 3 h.

2.3 Simulation of the Humidifying Process

When a porous medium gradually humidifies from a completely dry state, the evolution of the distributions of pore water and air can be approximately described as follows: First, water forms liquid bridges over the narrow spaces between pores and fills some of the dead-end pores in the medium. Pore water forms films on the surfaces of the granules. As the degree of water saturation increases, the liquid bridges widen, and the thickness of the water films increases. Then, the liquid bridges gradually connect to each other. When the proportion of air is further reduced, the connections between the bubbles are blocked, and eventually only isolated bubbles exist in the pores.

Figure 2 shows a 2D-CT (computed tomography) scanning image as an example. As shown in Figure 2, the humidifying process is achieved approximately in accordance with the description. This process can be simply divided into two steps: first, the gas particles are sorted according to degree of influence from the solid surface; second, these gas particles are converted into water particles in a certain order. The first step adopts the particle pairing algorithm of the SPH method. For the dry state, there are both solid particles and gas particles in the identification zone of a gas particle located near the surface of the granules (Figure 3). The parameter p is defined as the proportion of solid particles in the identification zone. In Figure 3, the p value for point A is greater than that for point B, and the p value for point B is greater than that for point C. According to their evolution, these particles located in dead-end pores (such as particle A) or near the surface of the solid phase (such as particle B) are wetted before those located in the centers of the pore channels (such as particle C). Moreover, the p values of the former are higher than those of the latter. Hence, the wetting order of the gas particles is basically in agreement with the order of the p values.

FIGURE 2. Simulation of the humidifying process for a 2D porous medium obtained by processing a 2D-CT scanning image (The blue particles denote water particles). (A) S_r =0.24 (B) S_r =0.53 (C) Sr =0.83.

FIGURE 3. The identification zone is similar to the influence domain of a particle, and the radius is defined as a multiple of the smoothing length.

Similarly, the humidifying process of SAP1 can also be obtained by the algorithm, as shown in Figure 4.

FIGURE 4. Simulation of the humidifying process for a 3D granular medium (Note: the corner is removed to reveal the internal details). (A) S_r =0.12 (B) S_r =0.40 (C) Sr =0.84.

2.4 Simulation Scheme

The steady-state method is a conventional laboratory measurement method for modeling the thermal conductivity. When the system reaches equilibrium, the heat transfer rate is equal to the heat dissipation rate in a stable thermal conduction process. According to Fourier’s law, the heat flux q of an ideal infinite plate satisfies (Cundall and Strack, 1979; Bai et al., 2019).

where ΔT is the temperature difference between the two sides of the plate, and L is the thickness of the plate.

In this calculation, an infinite plate is simulated by applying periodic boundary conditions to a square or cubic medium. The boundary conditions for the 2D simulation experiment are applied as shown in Figure 5. The upper and lower edges of the square medium are set as periodic boundaries so that the square medium can be periodically extended into an infinite plate with a certain thickness L. The left and right boundaries of the medium are set at different temperature values acting as the initial boundary conditions. For SPH particle sets, if the temperature difference is T₀, the average heat flux is calculated as

FIGURE 5. Diagram showing the boundary conditions of the steady-state simulation experiment.

When the temperature field reaches a steady state, the thermal conductivity of the medium can be determined by Eq. 12 ($\lambda =\raisebox{1ex}{$\overline{q}L$}\!\left/ \!\raisebox{-1ex}{$T_{\text{0}}$}\right.$).

The application of the finite difference method in SPH has been proposed in Bai’s work (Bai et al., 2018). In this work, the differential equations are applied in the post processing step of the SPH simulation experiments to output the temperature gradient. For some complex fluid-structure coupling boundary problems, some scholars (Huang et al., 2016; Liu and Zhang, 2019; Ye et al., 2019) have also proposed many improved calculation methods that can further improve the calculation accuracy (Liu and Liu, 2010).

3 Verification Using Typical Cases

3.1 Case 1: Homogeneous Medium

For this case, we assume that the infinite plate is a homogeneous medium with density ρ = 1000 kg/m³, specific heat capacity c = 1000 J/(kg°C), and default thermal conductivity λ₀ = 1 W/(m°C) (Bai et al., 2020; Bai et al., 2021a). The initial temperature of the particles is set to 0. In this calculation, the boundary conditions are applied as shown in Figure 5. The thickness of the square medium is L = 0.012 m. The other calculation parameters are listed in Table 1.

TABLE 1. Calculation parameters for the 2D case.

Figure 6A shows that the temperature distribution in the medium gradually becomes more linear. The dimensionless time factor is defined as t_d = κt/L² (Bai et al., 2020). Figure 6B shows that the calculated values of the thermal conductivity rapidly approach the baseline value during the dimensionless time t_d.

FIGURE 6. Calculation results of the homogeneous medium. (A) The temperature distribution (B) $\lambda$ -td trend line

3.2 Case 2: Three-Phase Layered Medium

As shown in Figure 7A, it is assumed that the solid, liquid, and gas components are arranged in layers, and each component occupies the same volume in the medium. The physical parameters of the materials are listed in Table 2. If boundaries 1 and 2 are treated as periodic boundaries and boundaries 3 and 4 are treated as the temperature boundaries, the three layers are connected in series. In contrast, if boundaries 3 and 4 are treated as periodic boundaries while boundaries 1 and 2 are treated as the temperature boundaries, the three layers are connected in parallel. In this simulation, the parameters are the same as those in Case 1 (Table 1) and are in agreement with the parameters in the existing literature (Alrtimi et al., 2014; Yang, 2019).

FIGURE 7. Three-phase layered medium and the calculation results. (A) Sketch map of the three-phase layered medium (B) Trend line of $\lambda$ and t_d.

TABLE 2. Thermal parameters of the materials involved in the calculation.

According to Eq. 1, the theoretical value of effective thermal conductivity for the three-phase series medium is 0.0745 W/(m°C), while that for the three-phase parallel medium is 2.742 W/(m°C).

Figure 7B shows that for both connection modes, the predicted values of the effective thermal conductivity precisely converge to the theoretical values over the chosen calculation time. The dimensionless time factor for the three-phase medium is defined as t_d = κt/L², with κ being the thermal diffusion coefficient of water (Bai et al., 2021c). The results show that the proposed simulation experiment is suitable for predicting the effective thermal conductivities of multiphase media.

4 The Relationship Between the Normalized Thermal Conductivity and Saturation Degree

In the SAP1 calculation, the temperature boundary is set to T₀ at the upper boundary (plate z = L) and set to 0 at the lower boundary (plate z = 0), while the other boundaries are set as periodic boundaries. The granular medium is simulated as an infinite porous layer with a thickness equal to L. The physical parameters for the components are assigned according to Table 2. The time step Δt is set to 0.0001 s. The simulated steady-state dimensionless temperature distribution of the sample in Figure 4 is predicted (Figure 8A).

FIGURE 8. 3D stacked media and calculation results. (A)Temperature distribution at a steady state (B) Predicted thermal conductivity over time.

Figure 8B shows that the effective thermal conductivity value gradually stabilizes over time. To ensure the precision of the three significant figures in the predicted λ value, the calculation process continues until the fourth significant digit stops changing in the time interval of Δt_d = 0.01. As shown in Figure 8B, the final result for the effective thermal conductivity for the medium reaches 2.75 W/(m°C), in agreement with the Wiener limits [the upper bound being 4.32 W/(m°C) and the lower bound being 0.06 W/(m°C), as calculated by Eq. 1].

A sample named SAP2 (Figure 9), with a granule radius of a = 0.8 mm, is generated as the control sample. SAP2 has the same porosity as SAP1. By applying the humidifying algorithm to both SAP1 and SAP2, two sets of samples with incremental degrees of saturation are obtained. In the simulation, the solid component for each group of soil samples is assigned first as quartz and then as feldspar. Thus, four sets of simulation experimental results are obtained. For feldspar, the values of λ, ρ and c are 2.3 W/(m°C), 2650 kg m⁻³ and 710 J/(kg°C), respectively. The trend line for λ_r versus S_r is plotted in Figure 10.

FIGURE 9. Particle distribution. (A) The control sample (named SAP2); (B) the discretization of SAP2 in the dry state.

FIGURE 10. The relationship between λ_r and S_r. (A) Trend line of λ_r versus S_r (B) Linear relation between 1/λ_r and 1/S_r.

The results show that the slope of the λ_r-S_r curve reaches a maximum when the degree of saturation is 0 and then decreases gradually. This curve indicates that the formation of a liquid bridge in the early stage of the humidifying process can improve the thermal conductivity of porous media. This phenomenon is consistent with the conclusions from other reference (Bai et al., 2019). Figure 10B shows that 1/λ_r has a good linear relationship with 1/S_r, which agrees with the empirical rule in Eq. 4. In addition, a comparison of the solid and hollow data points in Figure 10 shows that the slope of the linear fit for the 1/λ_r−1/S_r graph varies for different materials and granule sizes. Thus, the influence coefficient α from Eq. 4 is related to the thermal conductivity value of the solid phase and the granularity of the medium.

5 Conclusion

This research deduces a pore-scale simulation method to predict the effective thermal conductivities of unsaturated porous media using the SPH method. This study also discusses the relationship between the thermal conductivity and the degree of saturation in porous media. A concise method to simulate the wetting process of porous media is proposed based on pore-scale SPH simulation.

Based on a dry medium, the humidifying process can be modeled by sorting the gas particles according to the influence degree of the solid surface and then converting these gas particles into water particles in a certain order. This approach models the formation of liquid bridges and water films on granules. For three-phase layered media, the thermal conductivity values predicted by the simulation experiments agree with the theoretical values. This shows that the simulation experimental method is accurate for multiphase media. The particle fluid code is employed to generate the granular samples on which the simulation experiments are implemented.

The simulation experiments enable an investigation into the laws governing the relationship between the thermal conductivity and degree of saturation. The results show that the reciprocal of the normalized thermal conductivity and the reciprocal of the degree of saturation are linearly related. This result is consistent with the empirical formula proposed by Cote and Konrad. This verification of the law provides another way to demonstrate the usefulness and credibility of the humidifying process.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

The main contribution of BB in this paper is methodology, and the other authors (YW, DR, and FB) contributed in investigation and analysis.

Funding

This research was funded by the Beijing Natural Science Foundation (822023) and National Natural Science Foundation of China (52079003).

Conflict of Interest

YW was employed by the company Beijing Uni.-Construction Group Co., Ltd. DR was employed by the company Guangzhou Expressway Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Abbreviations

a, granule radius; c, specific heat capacity; f, field function; h, smoothing length of the kernel function; i, serial number of objective particles; j, serial number of particle i’s adjacent particles; L, thickness of the medium; m, mass of a particle; N, total number of particles within the influence area of particle i; p, the proportion of the solid particles in the identification zone; r, particle distance; R, relative distance; Sr, degree of saturation; t, time; t_d, dimensionless time, t_d = D₀t/L²; T, temperature; x, x’, x_i, x_j, distance vector of particle; W, kernel function; α, influence coefficient; α_d, normalized constant in the kernel function; $\nabla$, vector gradient operator; κ, thermal diffusion coefficient, equal to λ/ρc; λ, thermal conductivity; λ_r, normalized thermal conductivity; Δx, particle spacing; ρ, density of particle; $\phi$, porosity; $\varphi$, volume proportion of each component; Ω, integral domain.

References