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: f ( x ) ≈ f ˜ ( x ) = ∫ Ω f ( x ′ ) W ( | x − x ′ | , h ) d x ′ , (6) 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): W ( R , h ) = α d × { 2 / 3 − R 2 + R 3 / 2   0 ≤ R < 1 ( 2 − R ) 3 / 6         1 ≤ R < 2 0                 R ≥ 2 (7) 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 f ( x i ) ≈ ∑ j = 1 N m j ρ j f ( x j ) W ( | x i − x j | , h ) , (8) 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). ∂ T ∂ t = λ ρ c ⋅ ∇ 2 T (9) 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). ( ∂ T ∂ t ) i = 1 ρ i c i ∑ j = 1 N m j ρ j ⋅ 4 λ i λ j   ( λ i + λ j )   ( T i − T j ) ⋅ x i − x j r i j 2 ∇ i W i j (10) where r _ij =| x _ij |.

The size of each time interval is suggested to be (Bai et al., 2018). Δ t ≤ min ( 0.15 h 2 κ ) (11) where κ is the thermal diffusion coefficient equal to λ/ρc.

The explicit Euler method is employed for the integration step, f ˜ ( t + Δ t ) = f ( t ) + Δ t d f d t | 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).

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

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.

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.

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

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). q = λ Δ T L (12) 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 q ¯ = 1 Ν ∑ i = 1 Ν λ i ( d T d x ) i (13)

When the temperature field reaches a steady state, the thermal conductivity of the medium can be determined by Eq. 12 ( λ = q ¯ L T 0 ).

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

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.

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.

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

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.