Problem 4 - 3D Thermal Consolidation Test, Coupling with FLAC3D Simulator (by Guanlong Guo)

This example is contributed by Dr. Guanlong Guo, Energy Geosciences Division , LBNL , using TOUGH4 build No.: 2407301S

TOUGH4-FLAC3D is a two-way sequential coupling simulator, as detailed in Section 4. The coupling between the codes is achieved by iteratively writing/reading the binary files generated by each code during the simulation process. At the beginning of each time step, the porosity is updated based on the increments of volumetric strain, pore pressure and temperature in the previous time step. Then, the updated porosity is further corrected by the increment of the pressure and the temperature during the iteration process at the current time step. Mathematically, the porosity is updated as (Kim et al., 2012; Tounsi et al., 2023)

$\phi^{i+1}-\phi^i=\frac{{{(B}^i)}^2+\left(B^i-\phi^i\right)\left(1-B^i\right)}{K_{dr}^i}\sum_{\alpha=1,2}\left[S_\alpha\left(p_\alpha^{i+1}-p_\alpha^i\right)\right]+\phi^i\alpha_T\left(T^{i+1}-T^i\right)-\Delta\phi_c$ (10-6)

$\Delta\phi_c=\frac{\left(B^i\right)^2}{K_{dr}^i}\sum_{\alpha=1,2}\left[S_\alpha\left(p_\alpha^i-p_\alpha^{i-1}\right)\right]+B^i\alpha_T\left(T^i-T^{i-1}\right)-B^i\left(\epsilon_v^i-\epsilon_v^{i-1}\right)$ (10-7)

where B is the Biot's coefficient, $K_{dr}$ is the drained bulk modulus, $S_\alpha$ is the degree of saturation of phase $\alpha$ , $\alpha_T$ is the volumetric thermal expansion coefficient, T is the temperature and $\epsilon_v$ is the volumetric strain.

Once the coupled TH process reaches equilibrium, the equivalent pore pressure is calculated and then transferred to FLAC3D. The coupling discussed above is the direct coupling between the TH process and the mechanical process. In addition, some indirect couplings, such as changes in permeability and capillary pressure, have been implemented in the coupling code.

The thermal consolidation problem in this section assumes a point heat source buried in a saturated porous media. The analytical solution of the thermal consolidation problem is referenced in (Booker & Savvidou, 1985; Xu et al., 2020). Figure 10-6 presents the geometry and meshing strategy of the simulation domain. Due to the symmetric property, the domain below (50m × 50m× 50m) accounts for 1/8 of the entire domain. A heat source, 50 J/s, is applied to the element at the center of the whole domain, i.e. the element at the left-bottom corner of the domain in Figure 10-6. The initial isotropic stress is -0.1 MPa. The initial pore pressure is 0.1 MPa. The initial temperature is 25 °C. Roller boundaries are applied to the symmetric boundaries, while the other boundaries are subject to a compressive stress, -0.1 MPa. No flux boundary is applied to the symmetric boundaries and a constant pressure (0.1 MPa) is applied to the other boundaries.

Figure 10-6 X-Z cross-section of the 3D simulation domain

The material properties for the 3D thermal consolidation are given in Table 10-1.

Table 10-1 Material properties for the 3D thermal consolidation test

Figure 10-7 compares the analytical solution with the simulations results at a monitoring point M (1.2m,1.2m,1.2m). The simulation results are in good agreement with the analytical solution for temperature, pore pressure, shear stress. Slight deviations are observed in the simulation results for normal stress and displacement, as seen in Figure 10-7 (c) and (d). This difference may be attributed to the transfer of different time-stepping strategies in the simulation for enhancing computational efficiency. Using small time steps throughout the whole simulation may help prevent these discrepancy issues. In general, the TOUGH4-FLAC coupling code can accurately simulate the coupled THM problems.

Figure 10-7 Comparison between simulated results with analytical solutions

Input files

Output files