Hybrid fracture model deep fracture reservoir thermal flow coupling numerical simulation method and system

Abstract

The application discloses a kind of deep fracture reservoir thermal flow coupling numerical simulation method and system based on mixed fracture model, and relates to fractured reservoir numerical simulation technical field.The method constructs three types of grid systems of large-scale fracture, natural fracture and matrix in the same calculation domain, respectively using embedded discrete fracture model and double-pore double-permeability model to be coupled representation.Through establishing four types of grid connection relationship, flow and thermal conductivity are calculated differentially, especially introducing improved Vermeulen function to carry out non-steady-state dynamic correction to the exchange between natural fracture and matrix.Based on finite volume method, fully implicit residual equation set is established and solved, and multi-field simulation results are output.The application considers the representation accuracy and calculation efficiency of multi-scale medium under the premise of not multi-level profiling of matrix, and can provide technical support for deep oil and gas reservoir injection-production development and geothermal field efficient exploitation.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for fractured reservoirs, specifically to a method and system for thermal-fluid coupling numerical simulation of deep fractured reservoirs using a hybrid fracture model. Background Technology

[0002] Deep fractured reservoirs typically develop both natural and artificial fractures, and can form large-scale artificial fractures after fracturing. Therefore, they exhibit multi-scale media characteristics, with large-scale fractures, natural fractures, and matrix coexisting. During water injection development, cold fluid heating, and geothermal circulation projects, fluids preferentially migrate along high-conductivity fractures and continuously exchange mass and heat with natural fractures and the matrix. This results in significant heterogeneity and multi-scale coupling characteristics in the evolution of pressure, saturation, and temperature.

[0003] Among existing methods for simulating fractured reservoirs, embedded discrete fracture models can characterize fracture geometry without body-fitting the background mesh, making them suitable for handling complex, large-scale fractures. However, they typically do not separate the natural fracture system and the matrix system into a dual-pore, dual-permeability exchange system. Dual-pore, dual-permeability models can describe the macroscopic exchange between natural fractures and the matrix, but they usually cannot adequately describe large-scale fractures, and the natural fracture-matrix exchange often uses quasi-steady-state conductivity. Multi-medium models can improve the accuracy of matrix-natural fracture exchange characterization, but they require multi-level matrix subdivision, significantly increasing the number of meshes and computational costs.

[0004] On the other hand, while some heat-fluid coupling simulation methods can consider the coupling between the temperature field and the seepage field, they typically do not simultaneously distinguish between large-scale cracks, natural cracks, and the matrix at the macroscopic engineering scale, nor do they construct flow conductivity and thermal conductivity separately for the different connectivity relationships between the three types of media. Therefore, directly superimposing existing heat-fluid coupling equations onto conventional crack models makes it difficult to simultaneously achieve multi-scale characterization accuracy and engineering computational efficiency.

[0005] Therefore, there is an urgent need for a numerical simulation method for thermal-fluid coupling of deep fractured reservoirs that can explicitly characterize large-scale fractures, depict the dual-pore and dual-permeability behavior of natural fractures and matrix under the same computational framework, and simultaneously calculate flow conductivity and thermal conductivity for different grid connection relationships. In particular, it is necessary to perform conductivity correction for early unsteady mass exchange and heat exchange between natural fractures and matrix. Summary of the Invention

[0006] This invention provides a numerical simulation method and system for thermal-fluid coupling in deep fractured reservoirs using a hybrid fracture model. It solves the technical problem of simultaneously and accurately characterizing the fluid and heat exchange of large-scale fractures, natural fractures, and the matrix at a macroscopic engineering scale, while also addressing the inability to balance simulation accuracy and computational efficiency.

[0007] This invention is achieved through the following technical solution:

[0008] Firstly, this application provides a numerical simulation method for thermal-fluid coupling in deep fractured reservoirs using a hybrid fracture model, comprising:

[0009] Obtain basic parameters of the study area;

[0010] Based on the aforementioned fundamental parameters, a computational domain is established, and within the same computational domain, a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system are constructed. The large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, while the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model.

[0011] Multiphase flow mass conservation equations and energy conservation equations are established in the matrix grid system, natural fracture grid system, and large-scale fracture grid system, respectively.

[0012] Four types of mesh connection pairs are established, and the flow conductivity and thermal conductivity are calculated for each type of mesh connection pair. The conductivity of the fourth type of connection pair is dynamically corrected based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity.

[0013] Based on the corrected flow conductivity and thermal conductivity, the mass conservation equation and energy conservation equation of the multiphase flow are discretized by finite volume to establish a fully implicit residual equation system.

[0014] The Newton-Raphson iterative method was used to solve the fully implicit residual equations, and the numerical simulation results of thermal-fluid coupling for each grid system at different time steps were output.

[0015] A further optimization scheme is that the fluid parameters in the basic parameters include fluid density, viscosity, compressibility coefficient, relative permeability curve, specific heat capacity, and thermal conductivity coefficient;

[0016] The reservoir geological parameters include reservoir depth, thickness, initial pressure, and initial temperature;

[0017] The rock physical properties include porosity, permeability, rock density, rock specific heat capacity, and rock thermal conductivity.

[0018] The parameters of the natural fractures include the natural fracture spacing, porosity, permeability, specific heat capacity, and thermal conductivity.

[0019] The parameters of the large-scale cracks include crack geometric coordinates, crack length, crack height, crack width, permeability, porosity, specific heat capacity, and thermal conductivity.

[0020] A further optimization scheme is that the matrix grid system and the natural fracture grid system have a spatial correspondence, and the two are allowed to have fluid flow and heat transfer between adjacent grids within their respective systems through a dual-pore dual-permeability model.

[0021] A further optimization scheme is proposed for the four types of mesh connection pairs:

[0022] The flow conductivity and thermal conductivity of the first type of connection pair are calculated based on the interface area of ​​the connected grids, the distance from the grid center to the grid interface, the absolute permeability, and the effective thermal conductivity.

[0023] The flow conductivity and thermal conductivity of the second type of connection pair are calculated based on the area obtained by the large-scale fracture grid being cut by the matrix / natural fracture grid, the width of the large-scale fracture, the average vertical distance from the natural fracture / matrix grid to the large-scale fracture grid, the absolute permeability, and the effective thermal conductivity.

[0024] The flow conductivity and thermal conductivity of the third type of connection pair are calculated based on the number of intersecting large-scale crack meshes, the number of edges of the large-scale crack meshes, the edge length of the large-scale crack meshes, the distance from the large-scale crack meshes to the edges, the absolute permeability, and the effective thermal conductivity.

[0025] The fourth type of connection relationship is based on the improved Vermeulen function to calculate the flow conductivity correction coefficient and thermal conductivity correction coefficient as a function of pressure or temperature, and the corrected natural crack-matrix flow conductivity and thermal conductivity are obtained respectively.

[0026] A further optimization scheme is to first calculate the flow conductivity and thermal conductivity based on the mesh volume, natural crack spacing, absolute permeability, and effective thermal conductivity for the fourth type of connection pair.

[0027] Then, based on the improved Vermeulen function, the flow conductivity correction coefficient and thermal conductivity correction coefficient for pressure and temperature changes are calculated, and the flow conductivity and thermal conductivity are corrected respectively with the correction coefficients to obtain the corrected flow conductivity and thermal conductivity of the fourth type of connection pair.

[0028] A further optimization involves using an improved Vermeulen function to correct the conductance of the fourth type of connection pairs. The corrected conductance of the fourth type of connection pairs is:

[0029] ;

[0030] In the formula, TI ij,final γi represents the corrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid; γ1 and γ2 are empirical coefficients used to fit unsteady mass or heat transfer processes; δi iδi represents the pressure or temperature of the i-th matrix mesh; δ0 represents the initial pressure or temperature of the mesh at the start of the simulation; δi ... j The pressure or temperature of the j-th natural fracture grid; TI ij The uncorrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid;

[0031] A further optimized scheme is that the fully implicit residual equation set includes phase quality residual equations and energy residual equations, specifically:

[0032] ;

[0033] Where, when χ=β, it represents the β-phase quality residual equation:

[0034] ;

[0035] Where, when χ=T, it represents the energy residual equation:

[0036] ;

[0037] In the formula, R χi ρ represents the mass or energy residual in the i-th grid. βi,j λ is the density of phase β at the interface between the i-th and j-th grids; βi,j TI represents the mobility (permeability / viscosity) of the β phase at the interface between the i-th and j-th grids. ij p represents the flow or conductivity between connected grids. W W represents the bottom hole pressure. i V is the well index of the i-th grid; i Let be the volume of the i-th grid; Δt is the time step.

[0038] A further optimization scheme is that the multiphase flow mass conservation equation includes the flow term, accumulation term, and well source / sink term for each phase fluid in the corresponding grid system; the energy conservation equation includes the convection heat transfer term, heat conduction term, heat storage term, and well source / sink heat term.

[0039] Secondly, this application provides a hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation system, comprising:

[0040] The basic parameter acquisition module is used to acquire basic parameters of the study area;

[0041] The mesh construction module is used to establish a computational domain based on the basic parameters, and to construct a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system within the same computational domain; wherein, the large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, and the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model;

[0042] The model building module is used to establish the multiphase flow mass conservation equation and energy conservation equation in the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system, respectively.

[0043] The connectivity calculation module is used to establish four types of mesh connectivity pairs, calculate the flow conductivity and thermal conductivity for each type of mesh connectivity pair, and dynamically correct the conductivity of the fourth type of connectivity pair based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity.

[0044] The residual equation system establishment module is used to perform finite volume discretization on the multiphase flow mass conservation equation and energy conservation equation based on the corrected flow conductivity and thermal conductivity, and establish a fully implicit residual equation system.

[0045] The iterative solution module is used to solve the fully implicit residual equations using the Newton-Raphson iterative method and output the thermal-fluid coupling numerical simulation results of each grid system at different time steps.

[0046] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0047] By constructing a hybrid fracture model that couples an embedded discrete fracture model with a dual-pore dual-permeability model, three types of media—large-scale fractures, natural fractures, and matrix—are collaboratively characterized within the same computational domain. Four types of mesh connections are established to calculate flow and thermal conductivity differently. In particular, an improved Vermeulen function is introduced to perform unsteady-state dynamic correction for the exchange between natural fractures and the matrix. This approach avoids multi-level redundant subdivision of the matrix while accurately depicting the fluid transport and heat transfer patterns within multi-scale media. It significantly reduces computational costs and improves the prediction accuracy of pressure, saturation, and temperature fields in the simulation of injection, production, and geothermal circulation in deep fractured reservoirs. This model is applicable to efficient oilfield exploitation technologies. Attached Figure Description

[0048] To more clearly illustrate the technical solutions of the exemplary embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0049] In the attached image:

[0050] Figure 1 This is a flowchart illustrating the implementation of the hybrid fracture model deep fracture reservoir thermal flux coupling numerical simulation method provided in this application embodiment;

[0051] Figure 2 The matrix mesh pressure variation curve obtained from the simulation calculation of Example 1 provided in this application;

[0052] Figure 3 The temperature change curve of the matrix mesh obtained by simulation calculation in Example 1 provided in this application;

[0053] Figure 4 A schematic diagram of the physical model of Embodiment 2 provided in this application;

[0054] Figure 5 The oil saturation distribution map obtained by simulation calculation in Example 2 provided in this application;

[0055] Figure 6 The temperature distribution map obtained by simulation calculation for Embodiment 2 provided in this application;

[0056] Figure 7 A bar chart comparing the number of meshes and computation time for different models in Embodiment 2 provided in this application;

[0057] Figure 8 A schematic diagram of the physical model of Embodiment 3 provided in this application;

[0058] Figure 9 The temperature distribution map obtained by simulation calculation in Example 3 provided in this application;

[0059] Figure 10 A bar chart comparing the number of meshes and computation time for different models in Example 3 provided in this application. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0061] First, some of the technical terms used in this application will be explained to help those skilled in the art understand this application.

[0062] EDFM: Embedded Discrete Fracture Model;

[0063] DPDK: Dual Porosity Dual Permeability, a dual-pore, dual-permeability model.

[0064] Firstly, this application provides a numerical simulation method for thermal-fluid coupling in deep fractured reservoirs using a hybrid fracture model, comprising:

[0065] Step S1: Obtain the basic parameters of the study area;

[0066] Step S2: Establish a computational domain based on the aforementioned basic parameters, and construct a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system within the same computational domain; wherein, the large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, and the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model;

[0067] Step S3: Establish the multiphase flow mass conservation equation and energy conservation equation in the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system, respectively;

[0068] Step S4: Establish four types of mesh connection pairs, calculate the flow conductivity and thermal conductivity for each type of mesh connection pair, and dynamically correct the conductivity of the fourth type of connection pair based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity.

[0069] Step S5: Based on the corrected flow conductivity and thermal conductivity, the multiphase flow mass conservation equation and energy conservation equation are discretized by finite volume to establish a fully implicit residual equation set.

[0070] Step S6: Solve the fully implicit residual equations using the Newton-Raphson iterative method and output the numerical simulation results of thermal-fluid coupling for each grid system at different time steps.

[0071] This embodiment achieves collaborative characterization of three types of media—large-scale fractures, natural fractures, and matrix—within the same computational domain by coupling an embedded discrete fracture model with a dual-pore, dual-permeability model. By establishing four types of mesh connectivity and calculating conductivity differently, and particularly by utilizing an improved Vermeulen function to perform unsteady dynamic correction on the exchange between natural fractures and the matrix, it balances multi-scale characterization accuracy and computational efficiency without performing multi-level redundant matrix subdivision. This significantly improves the accuracy of deep reservoir thermal-fluid coupling simulation and is suitable for efficient oilfield exploitation technologies.

[0072] In one embodiment, the basic parameters include fluid parameters, reservoir geological parameters, rock physical property parameters, natural fracture parameters, and large-scale fracture parameters; Step S1: Obtain the basic parameters of the study area, specifically including the following steps:

[0073] Step S11: Collect fluid property information of the study area and integrate it to obtain a set of fluid parameters; specifically, the fluid parameters include fluid density, viscosity, compressibility coefficient, relative permeability curve, specific heat capacity and thermal conductivity.

[0074] Step S12: Extract geological background data of the study area and compile a set of reservoir geological parameters; specifically, the reservoir geological parameters include reservoir depth, thickness, initial pressure and initial temperature.

[0075] Step S13: Determine the physical properties of the rocks in the study area and compile a set of rock physical property parameters; specifically, the rock physical property parameters include porosity, permeability, rock density, rock specific heat capacity and rock thermal conductivity.

[0076] Step S14: Statistically analyze the development characteristics of natural cracks in the study area and summarize the set of natural crack parameters; specifically, the natural crack parameters include natural crack spacing, porosity, permeability, specific heat capacity and thermal conductivity.

[0077] Step S15: Collect large-scale fracture geometry and physical property data in the study area to construct a set of large-scale fracture parameters; specifically, the large-scale fracture parameters include fracture geometric coordinates, fracture length, fracture height, fracture width, permeability, porosity, specific heat capacity, and thermal conductivity.

[0078] This embodiment collects relevant parameters of fluid, geology, rock, natural fractures and large-scale fractures in the study area by hierarchical classification, forming a complete set of basic parameters. This provides accurate input basis for subsequent computational domain construction, grid system division and establishment of multiphase flow and energy conservation equations, ensuring the reliability and prediction accuracy of multi-scale heat-fluid coupling simulation.

[0079] In one embodiment, step S2: A computational domain is established based on the fundamental parameters, and a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system are constructed within the same computational domain; wherein, the large-scale fracture mesh system is explicitly characterized using an embedded discrete fracture model, and the natural fracture mesh system and the matrix mesh system are characterized using a dual-pore dual-permeability model, specifically including the following steps:

[0080] Step S21: Delineate the simulation space range based on the acquired basic parameters of the study area and generate a computational domain covering the target reservoir.

[0081] Specifically, based on the burial depth, thickness, and planar distribution range of the reservoir geological parameters, the three-dimensional spatial boundary of the computational domain is determined to ensure that it includes all matrix, natural fractures, and large-scale fracture areas to be simulated.

[0082] Step S22: Simultaneously construct a matrix mesh system and a natural crack mesh system that are spatially completely overlapping within the computational domain.

[0083] Specifically, the matrix region is uniformly or non-uniformly divided to generate a matrix mesh system; simultaneously, a natural crack mesh system with identical geometric dimensions and quantity is generated at the same spatial location.

[0084] Among them, the matrix grid system represents matrix rock blocks with low porosity and low permeability, which mainly serve as fluid storage space; the natural fracture grid system represents high-porosity and high-permeability high-speed seepage channels; the two achieve functional separation of "matrix storage and fracture conduction" by endowing them with completely different rock properties, and flow and heat transfer can occur between adjacent grids within both systems.

[0085] Step S23: Explicitly construct a large-scale crack mesh system within the computational domain using an embedded discrete crack model.

[0086] Specifically, based on the geometric coordinates, length, height, and width of the large-scale crack parameters, the large-scale cracks are cut into the background matrix mesh using a non-body-fitting method, forming a large-scale crack mesh system that coexists with the matrix mesh system and the natural crack mesh system. The geometry of this mesh system is consistent with the actual cracks, representing the main channels with high conductivity and directly controlling the macroscopic migration path of the fluid.

[0087] This embodiment constructs three types of grid systems collaboratively within the same computational domain, which not only preserves the explicit geometric features of large-scale cracks but also achieves dual-pore dual-permeability coupling between natural cracks and the matrix. This avoids the multi-level partitioning of the matrix by the multi-medium model, providing an accurate grid foundation for the subsequent establishment of multiphase flow and energy conservation equations, and ensuring efficient characterization of flow and heat transfer between multi-scale media.

[0088] In one embodiment, step S3: establishing the multiphase flow mass conservation equation and energy conservation equation in the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system respectively, specifically includes the following steps:

[0089] Step S31: Based on the medium characteristics of the three types of grid systems, construct the basic framework of the mass conservation equation for multiphase flow;

[0090] Specifically, for the pore structure and fluid occurrence state of the matrix, natural fractures and large-scale fractures, the flow terms, accumulation terms and well source and sink terms that need to be included in the equations are clarified to ensure that the equations can reflect the mass change law of multiphase fluids in different media.

[0091] Step S32: Based on the energy transfer characteristics of the three types of grid systems, construct the basic framework of the energy conservation equation;

[0092] Specifically, by combining the differences in heat capacity between fluids and rocks and the heat transfer methods, the convective heat transfer term, heat conduction term, heat storage term, and well source heat sink term that need to be included in the equation are determined to ensure that the equation can reflect the energy accumulation and dissipation process in different media.

[0093] Step S33: Based on the basic framework of the constructed multiphase flow mass conservation equation, Darcy's law is used to describe the flow of fluid in porous media. Combining phase saturation constraints and the principle of mass balance, the specific expression of the multiphase flow mass conservation equation is derived, as shown in the following equation:

[0094]

[0095] In the formula, ρ β k is the density of the β phase; k is the absolute permeability; k rβ The relative permeability of the β phase; μ β The viscosity of the phase; p β q represents the pressure of the β phase; g represents gravitational acceleration; D represents the elevation; β for Phase well flow rate; () W Represents the source and sink term; ϕ is porosity; s β The saturation of the β phase.

[0096] Step S34: Based on the basic framework of the constructed energy conservation equation, and comprehensively considering the internal energy carried by fluid flow, the thermal energy stored in the rock skeleton, and the heat conduction driven by the temperature gradient, the specific expression of the energy conservation equation is derived, as shown in the following equation:

[0097]

[0098] In the formula, h β Enthalpy of the β phase; κ t ρ is the effective thermal conductivity; T is the temperature; (ρC) p ) t This is the effective volumetric heat capacity.

[0099] Step S35: Construct a parameter calculation system that includes the energy conservation equations for enthalpy, effective thermal conductivity, and effective volumetric heat capacity.

[0100] Enthalpy h:

[0101]

[0102] In the formula, h is the enthalpy; C p ρ is the specific heat capacity at constant pressure; T is the temperature; p is the pressure.

[0103] Effective thermal conductivity κt:

[0104]

[0105] In the formula, κ r κ is the thermal conductivity of the rock. βis the thermal conductivity of the β phase.

[0106] Effective volumetric heat capacity (ρC) p ) t :

[0107]

[0108] In the formula, ρ r C is the density of the rock. pr C is the specific heat capacity of the rock. pβ is the specific heat capacity of the β phase.

[0109] This embodiment establishes the basic framework of mass and energy conservation through steps S31 and S32, and derives the multiphase flow mass conservation equation by introducing Darcy's law in step S33. It then derives the energy conservation equation by integrating heat conduction and convection mechanisms in step S34, thus comprehensively describing the transport and heat transfer laws of fluids in the matrix, natural fractures, and large-scale fractures. Step S35 clarifies the calculation methods of key auxiliary parameters such as enthalpy, effective thermal conductivity, and effective volumetric heat capacity, providing necessary physical property parameter support for the aforementioned equation set. Based on the above steps, a multiphysics coupled mathematical model is constructed, realizing efficient simulation of multiphase fluids in complex fractured reservoirs.

[0110] In one embodiment, step S4: Establish four types of mesh connection pairs, calculate the flow conductivity and thermal conductivity for each type of mesh connection pair, and dynamically correct the conductivity of the fourth type of connection pair based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity. This specifically includes the following steps:

[0111] Step S41: Based on the constructed matrix grid system, natural fracture grid system, and large-scale fracture grid system, clarify the connection relationship types between different media and generate a list of four types of grid connection relationship pairs;

[0112] Specifically, the four types of connection pairs include the first type of connection pair consisting of adjacent grids within the same grid system, the second type of connection pair consisting of matrix grids or natural fracture grids and their internal large-scale fracture grids, the third type of connection pair consisting of intersecting large-scale fracture grids within the same grid system, and the fourth type consisting of natural fracture grids and matrix grids occupying the same spatial position and having completely overlapping geometric shapes within the framework of the dual-pore dual-permeability model.

[0113] Due to the significant difference in permeability between the natural fracture mesh and the matrix mesh, the fluid flow and heat transfer between them exhibit strong unsteady characteristics in the early stages. Traditional static conductivity cannot accurately describe this exchange process. Therefore, the fourth type of connectivity needs to be dynamically corrected based on the improved Vermeulen function to accurately characterize the fluid supply from the matrix to the fracture and the heat exchange mechanism between them.

[0114] Step S42: For each type of connection pair, calculate the flow conductivity and thermal conductivity by combining the geometric parameters and physical property parameters of the corresponding mesh.

[0115] Specifically, the general formulas for calculating the flow conductivity and thermal conductivity of the first type of connection pair are as follows:

[0116]

[0117] In the formula, TI ij For the flow or thermal conductivity of phase-connected matrix or natural crack mesh i and j; TI i TI represents the flow or thermal conductivity of the i-th grid. j Let be the flow or thermal conductivity of the j-th grid.

[0118] The general formulas for calculating the flow or thermal conductivity of the i-th and j-th grids are as follows:

[0119]

[0120] In the formula, ξ represents the absolute permeability or effective thermal conductivity of the i-th and j-th grids; A ij d is the area of ​​the interface between the i-th and j-th grids; i d is the distance from the center of the i-th grid to the interface; j Let be the distance from the center of the j-th grid to the interface.

[0121] Specifically, the general formula for calculating the flow conductivity and thermal conductivity of the second type of connection pair also adopts Equation (6), where the conductivity of the matrix / natural fracture mesh i and the large-scale fracture mesh j are respectively:

[0122]

[0123] In the formula, ξ is the absolute permeability or effective thermal conductivity of the i-th matrix / natural fracture grid; A Fj ω represents the area obtained by cutting the j-th large-scale fracture mesh by the i-th matrix / natural fracture mesh; Fj d represents the width of the j-th large-scale crack mesh; FjThis represents the average vertical distance from the i-th matrix / natural fracture grid to the j-th large-scale fracture grid.

[0124] Specifically, the general formulas for calculating the flow conductivity and thermal conductivity of the third type of connection pair are as follows:

[0125]

[0126] In the formula, n represents the number of intersecting large-scale crack meshes; z k Represents the number of edges in the k-th large-scale crack mesh; l km d represents the length of the m-th edge of the k-th large-scale crack mesh; km This represents the distance from the k-th large-scale crack mesh to the m-th edge.

[0127] Specifically, the general formulas for calculating the flow conductivity and thermal conductivity of the fourth type of connection pair are as follows:

[0128]

[0129] In the formula, V i Lx represents the mesh volume; Lx, Ly, and Lz represent the natural crack spacing.

[0130] Furthermore, by modifying the Vermeulen function to correct the flow conductivity and thermal conductivity of the fourth type of connection pairs, the general formulas for calculating the corrected flow conductivity and thermal conductivity are shown below:

[0131]

[0132] In the formula, TI ij,final γi represents the corrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid; γ1 and γ2 are empirical coefficients used to fit unsteady mass or heat transfer processes; δi i δi represents the pressure or temperature of the i-th matrix mesh; δ0 represents the initial pressure or temperature of the mesh at the start of the simulation; δi ... j The pressure or temperature of the j-th natural fracture grid; TI ij The uncorrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid;

[0133] This embodiment establishes four types of mesh connectivity and calculates flow and thermal conductivity differently. This invention achieves accurate quantification of the complex exchange behavior among three types of media: large-scale cracks, natural cracks, and the matrix.

[0134] The first type of connection uses the harmonic average method to calculate the conductance of adjacent grids, ensuring numerical stability within the computational domain;

[0135] The second and third types of connection pairs define exclusive conductivity algorithms for the embedded contact between matrix / natural cracks and large-scale cracks, as well as the cross contact between large-scale cracks, respectively, solving the problem of heat exchange calculation under complex crack geometry.

[0136] The fourth type of connectivity introduces an improved Vermeulen function to perform unsteady-state dynamic correction on the conductivity between natural cracks and the matrix, overcoming the limitation of traditional quasi-steady-state models that cannot accurately characterize the early rapid crossflow and heat transfer processes.

[0137] This design balances the accuracy of multi-scale medium characterization with computational efficiency without performing multi-level redundant partitioning of the matrix. It provides accurate physical input for the subsequent construction of fully implicit residual equations and significantly improves the predictive reliability of pressure, saturation, and temperature fields in deep fractured reservoirs.

[0138] In one embodiment, step S5: Based on the corrected flow conductivity and thermal conductivity, the multiphase flow mass conservation equation and energy conservation equation are discretized by finite volume to establish a fully implicit residual equation system, specifically including the following steps:

[0139] Specifically, the pressure, saturation, and temperature of each grid in the three types of grid systems are taken as unknowns, and the bottom hole flowing pressure, bottom hole temperature, and production rate are taken as solution variables. The fully implicit residual equation set includes the phase quality residual equation and the energy residual equation, as shown in the following equation.

[0140]

[0141] Where, when χ=β, it represents the β-phase quality residual equation:

[0142]

[0143] Where, when χ=T, it represents the energy residual equation:

[0144]

[0145] This embodiment clarifies four types of mesh connection relationships and calculates conductivity differently. It establishes a fully implicit residual equation system by combining finite volume discretization. This avoids the simplification of different medium connection relationships in traditional models. It also improves the accuracy of multi-scale thermal-fluid coupling simulation by modifying the Vermeulen function to correct the conductivity of natural cracks and matrix, while maintaining computational efficiency. This provides an accurate algebraic equation basis for subsequent solutions.

[0146] In one embodiment, step S6: Solving the fully implicit residual equations using the Newton-Raphson iterative method and outputting the thermal-fluid coupling numerical simulation results of each grid system at different time steps, specifically includes the following steps:

[0147] Step S61: Set the initial parameters for the solution process and prepare to enter the iterative calculation; specifically, the initial parameters include the time step, the maximum number of iterations, the residual norm convergence threshold, and the variable increment convergence threshold. These parameters are determined in advance according to the simulation requirements and computing resources.

[0148] Step S62: Based on the flow conductivity and thermal conductivity, the Newton-Raphson iterative solution is performed using the fully implicit residual equations. Specifically, the pressure, saturation, temperature, bottom hole flowing pressure, bottom hole temperature, and production rate at the current time step are taken as unknowns and substituted into the phase quality residual equation and energy residual equation to calculate the residual vector and Jacobian matrix. The unknowns are then updated using a linear solver. The phase quality residual equation and energy residual equation have been detailed in step S5.

[0149] Step S63: Determine whether the current iteration meets the convergence condition; specifically, calculate the norm of the residual vector and the absolute value of the increment of each unknown quantity. If both are less than the preset convergence threshold, the iteration is determined to be converged; otherwise, return to step S62 to continue iterating until the maximum number of iterations is reached.

[0150] Step S64: If the iteration converges, use the solution result of the current time step as the initial value of the next time step and update the time step to the next moment; specifically, assign the converged variables such as pressure, saturation, and temperature to the initial guess of the next time step, and the time step size can be adaptively adjusted according to the computational stability.

[0151] Step S65: Output the thermal-fluid coupling numerical simulation results of each grid system at the current time step; specifically, the results include the pressure field, saturation field, and temperature field of the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system, as well as dynamic parameters such as oil production, water production, gas production, bottom hole pressure, and bottom hole temperature of the production well.

[0152] This embodiment significantly improves the computational efficiency of thermal-fluid coupling simulation of deep fractured reservoirs by setting rigorous initial parameters and convergence criteria, and by adopting the Newton-Raphson iteration method combined with an adaptive time step strategy, while ensuring the numerical stability of the fully implicit residual equation system. In particular, by using the converged pressure, saturation, and temperature field results as the initial values ​​for the next time step, and dynamically outputting the multi-field distribution of three types of grid systems—matrix, natural fractures, and large-scale fractures—as well as the dynamic parameters of the production well, it achieves an effective correlation from micro-grid transmission to macro-development indicators, providing a high-precision and high-reliability quantitative basis for optimizing injection-production development schemes and predicting geothermal production capacity.

[0153] To verify the effectiveness of the hybrid fracture model deep fracture reservoir thermal flux coupling numerical simulation method provided by this invention, the following detailed explanation and analysis are provided in conjunction with specific implementation cases.

[0154] Example 1

[0155] Example 1 simulates the mass and heat transfer process of a single-phase (water-phase) micro-compressible fluid in a single three-dimensional matrix rock block. It only requires dividing the matrix into one matrix grid and one fracture grid, and treating the pressure and temperature of the fracture grid as constant external boundary conditions of the matrix grid. This verifies the rationality of the proposed correction calculation of flow conductivity and thermal conductivity between natural fractures and the matrix. Table 1 gives the parameters that can be used in Example 1.

[0156] Table 1. Basic Parameters Related to Example 1

[0157]

[0158] like Figure 2 and Figure 3 As shown, after the calculation is completed, the pressure and temperature changes of the matrix mesh can be output.

[0159] like Figure 2 As shown, the matrix pressure decreases rapidly in the early stages, eventually reaching the same value as the fracture pressure. Without considering corrections, the early decrease rate of the matrix pressure is relatively slow, resulting in a larger error compared to the accurate analytical solution in the early stages. Conversely, the calculation results obtained with the Vermeulen function correction and the improved Vermeulen function correction generally show high consistency with the accurate analytical solution. However, the improved Vermeulen function correction proposed in this invention provides higher accuracy than the original Vermeulen function correction. Figure 3 As shown, it can be observed that the rate of decrease in matrix temperature is slower than the rate of decrease in pressure. Without considering the correction, the matrix temperature change obtained by calculation differs significantly from the accurate result of the analytical solution. The Vermeulen function correction result has an error in the middle stage compared with the accurate result of the analytical solution. However, the improved Vermeulen function correction result proposed in this invention has a high consistency with the accurate result of the analytical solution.

[0160] Example 2

[0161] like Figure 4As shown, Example 2 uses a physical model with dimensions of 500×200×100 m. Three large-scale fractures are incorporated into the model, with a water injection well in the lower left corner and an oil production well in the upper right corner. The large-scale fractures are explicitly embedded using EDFM, and the natural fractures and matrix are characterized using DPDK. Flow conductivity and thermal conductivity are calculated simultaneously for all four types of connectivity relationships, according to the first type. Table 2 provides the parameters applicable to Example 2.

[0162] Table 2. Basic Parameters Related to Example 2

[0163]

[0164] like Figure 5 and Figure 6 As shown, after the calculation is completed, the oil saturation distribution and temperature distribution in the matrix grid system, the natural fracture grid system and the large-scale fracture grid system can be output.

[0165] like Figure 5 As shown, the distribution range of low oil saturation regions in the fracture system is significantly larger than that in the matrix system, indicating that injected water has stronger flow and displacement capabilities in the fracture system. This is because the overall permeability of the fracture system is higher, especially the large-scale fractures, which have strong conductivity, allowing injected water to preferentially migrate rapidly along the fracture system and form a large water-drive sweep area. Meanwhile, the oil saturation distribution on the right side shows a characteristic of extending along the direction of the large-scale fractures, indicating that the high-conductivity fractures have a significant control effect on the advancement direction of the water-drive front. In contrast, the matrix system has lower permeability, making it relatively difficult for injected water to enter and displace crude oil in the matrix pores. Therefore, the low oil saturation region in the matrix system is smaller, and the water-drive front mainly exhibits a non-uniform advancement characteristic controlled by large-scale fractures.

[0166] like Figure 6 As shown, the propagation range of the cryogenic front in the fracture system is significantly greater than that in the matrix system. This phenomenon is mainly due to the stronger seepage capacity of the fracture system, allowing injected cold water to migrate rapidly within the fracture network and propagate the cryogenic region further outward through convective heat transfer. Therefore, the cryogenic front in the fracture system extends a greater distance. In contrast, the matrix system has lower permeability, slower fluid transport speed, and its temperature changes are more influenced by conduction heat transfer between the fracture and matrix, as well as local fluid exchange. Consequently, the propagation range of the cryogenic front is relatively smaller.

[0167] like Figure 7As shown in Table 3, the calculation results of different fracture characterization models in Example 2 are compared. It can be seen that, compared with the traditional dual-pore dual-permeability model and the embedded discrete fracture model, the present invention increases the number of grids and the computation time, but makes up for the lack of explicit characterization of large-scale fractures, coupling of natural fractures and matrix with dual pores and dual permeability, and the ability to dynamically correct the fourth type of conductivity. Compared with the multi-medium model, it significantly reduces the number of grids and the computation time while maintaining the ability to characterize the multi-scale medium thermal-fluid coupling. In the cumulative oil production prediction of Example 2, the errors of the traditional dual-pore dual-permeability model, the multi-medium model, the embedded discrete fracture model, and the model of the present invention compared with the refined model are 10.98%, 6.88%, 7.91%, and 2.83%, respectively, with the model of the present invention having the smallest error.

[0168] Table 3 Comparison of Calculation Results of Different Crack Characterization Models in Example 2

[0169]

[0170] Example 3

[0171] Example 3 uses a physical model with dimensions of 310×310×100 m. The model includes two large-scale cracks, a water injection well in the lower left corner, and a heat extraction well in the upper right corner. Figure 8 As shown, large-scale cracks were explicitly embedded using EDFM, and natural cracks and the matrix were characterized using DPDK. Flow conductivity and thermal conductivity were calculated simultaneously for connections of types I to IV. Table 4 provides the parameters applicable to Example 3.

[0172] Table 4. Basic Parameters Related to Example 3

[0173]

[0174] After the calculation is completed, the temperature distribution in the matrix mesh system, the natural fracture mesh system, and the large-scale fracture mesh system can be output, such as... Figure 9 As shown.

[0175] like Figure 9 As shown, both the fracture system and the matrix system exhibit a temperature distribution characteristic that transitions from a low-temperature region at the injection end to a high-temperature region at the distal end. The low-temperature region in the fracture system is relatively large, and the overall temperature field shows a relatively continuous advancing characteristic, indicating that the injected cold water has a strong transport capacity and significant convective heat transfer within the fracture system. In contrast, the temperature distribution in the matrix system is influenced by both the large-scale fractures and the low permeability of the matrix. The advancing range of the low-temperature front is relatively small, and local temperature changes show some non-uniformity, indicating that heat transfer in the matrix system is mainly controlled by heat conduction between the fractures and the matrix, as well as local fluid exchange.

[0176] Table 5 lists a comparison of the calculation results of different crack characterization models in Example 3, such as... Figure 10 As shown, it can also be seen that compared with the traditional dual-pore dual-permeability model and the embedded discrete crack model, the present invention increases the number of grids and computation time, but makes up for the lack of explicit characterization of large-scale cracks, coupling of natural cracks and matrix in dual-pore dual-permeability, and the ability to dynamically correct the fourth type of conductivity. Compared with the multi-medium model, it significantly reduces the number of grids and computation time while maintaining the ability to characterize the multi-scale medium thermal-fluid coupling. In the output temperature prediction of Example 3, the errors of the traditional dual-pore dual-permeability model, the multi-medium model, the embedded discrete crack model, and the model of the present invention compared with the refined model are 19.91%, 9.00%, 12.45%, and 3.00%, respectively, and the model of the present invention has the smallest error.

[0177] Table 5 Comparison of Calculation Results of Different Crack Characterization Models in Example 3

[0178]

[0179] The above results demonstrate that the numerical simulation method for deep fractured reservoir thermal-fluid coupling based on a hybrid fracture model provided by this invention innovatively constructs a three-grid system with an embedded discrete fracture model coupled with a dual-pore dual-permeability model. It specifically uses the embedded discrete fracture model to explicitly characterize large-scale fractures, the dual-pore dual-permeability model to depict the natural fracture-matrix system, and dynamically corrects the flow conductivity and thermal conductivity of the fourth type of natural fracture-matrix connection relationship through an improved Vermeulen function. Example 1 shows that this correction can more accurately approximate the analytical solution and precisely capture the pressure and temperature changes in the matrix grid. Examples 2 and 3 show that the model maintains its multi-scale thermal-fluid coupling characterization capability while balancing simulation accuracy and computational efficiency.

[0180] Secondly, this application provides a hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation system, comprising:

[0181] The basic parameter acquisition module is used to acquire basic parameters of the study area;

[0182] The mesh construction module is used to establish a computational domain based on the basic parameters, and to construct a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system within the same computational domain; wherein, the large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, and the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model;

[0183] The model building module is used to establish the multiphase flow mass conservation equation and energy conservation equation in the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system, respectively.

[0184] The connectivity calculation module is used to establish four types of mesh connectivity pairs, calculate the flow conductivity and thermal conductivity for each type of mesh connectivity pair, and dynamically correct the conductivity of the fourth type of connectivity pair based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity.

[0185] The residual equation system establishment module is used to perform finite volume discretization on the multiphase flow mass conservation equation and energy conservation equation based on the corrected flow conductivity and thermal conductivity, and establish a fully implicit residual equation system.

[0186] The iterative solution module is used to solve the fully implicit residual equations using the Newton-Raphson iterative method and output the thermal-fluid coupling numerical simulation results of each grid system at different time steps.

[0187] The functions of each module in the above-mentioned hybrid fracture model deep fracture reservoir thermal flux coupling numerical simulation system correspond to the steps in the above-mentioned hybrid fracture model deep fracture reservoir thermal flux coupling numerical simulation method embodiment, and their functions and implementation processes will not be described in detail here.

[0188] Thirdly, embodiments of this application also provide a readable storage medium.

[0189] This application stores a hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation program on a readable storage medium, wherein when the hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation program is executed by a processor, it implements the steps of the hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation method as described above.

[0190] The method implemented when the hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation program is executed can refer to the various embodiments of the hybrid fracture model deep fracture reservoir thermal-fluid coupling numerical simulation method of this application, and will not be repeated here.

[0191] It should be noted that the sequence numbers of the embodiments in this application are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0192] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A numerical simulation method for thermal-fluid coupling in deep fractured reservoirs using a hybrid fracture model, characterized in that, include: Obtain basic parameters of the study area; Based on the aforementioned fundamental parameters, a computational domain is established, and within the same computational domain, a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system are constructed. The large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, while the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model. Multiphase flow mass conservation equations and energy conservation equations are established in the matrix grid system, natural fracture grid system, and large-scale fracture grid system, respectively. Four types of mesh connection pairs are established, and the flow conductivity and thermal conductivity are calculated for each type of mesh connection pair. The conductivity of the fourth type of connection pair is dynamically corrected based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity. Based on the corrected flow conductivity and thermal conductivity, the mass conservation equation and energy conservation equation of the multiphase flow are discretized by finite volume to establish a fully implicit residual equation system. The Newton-Raphson iterative method was used to solve the fully implicit residual equations, and the numerical simulation results of thermal-fluid coupling for each grid system at different time steps were output.

2. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 1, characterized in that, The fluid parameters in the basic parameters include fluid density, viscosity, compressibility, relative permeability curve, specific heat capacity, and thermal conductivity. The reservoir geological parameters include reservoir depth, thickness, initial pressure, and initial temperature; The rock physical properties include porosity, permeability, rock density, rock specific heat capacity, and rock thermal conductivity. The parameters of the natural fractures include the natural fracture spacing, porosity, permeability, specific heat capacity, and thermal conductivity. The parameters of the large-scale cracks include crack geometric coordinates, crack length, crack height, crack width, permeability, porosity, specific heat capacity, and thermal conductivity.

3. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 1, characterized in that, The matrix grid system and the natural fracture grid system are spatially corresponding, and the two systems are connected by a dual-pore dual-permeability model, which allows fluid flow and heat transfer between adjacent grids within their respective systems.

4. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 1, characterized in that, Among the four types of mesh connection relationships: The flow conductivity and thermal conductivity of the first type of connection pair are calculated based on the interface area of ​​the connected grids, the distance from the grid center to the grid interface, the absolute permeability, and the effective thermal conductivity. The flow conductivity and thermal conductivity of the second type of connection pair are calculated based on the area obtained by the large-scale fracture grid being cut by the matrix / natural fracture grid, the width of the large-scale fracture, the average vertical distance from the natural fracture / matrix grid to the large-scale fracture grid, the absolute permeability, and the effective thermal conductivity. The flow conductivity and thermal conductivity of the third type of connection pair are calculated based on the number of intersecting large-scale crack meshes, the number of edges of the large-scale crack meshes, the edge length of the large-scale crack meshes, the distance from the large-scale crack meshes to the edges, the absolute permeability, and the effective thermal conductivity. The fourth type of connection relationship is based on the improved Vermeulen function to calculate the flow conductivity correction coefficient and thermal conductivity correction coefficient as a function of pressure or temperature, and the corrected natural crack-matrix flow conductivity and thermal conductivity are obtained respectively.

5. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 4, characterized in that, For the fourth type of connection pair, the flow conductivity and thermal conductivity are first calculated based on the mesh volume, natural crack spacing, absolute permeability and effective thermal conductivity. Then, based on the improved Vermeulen function, the flow conductivity correction coefficient and thermal conductivity correction coefficient for pressure and temperature changes are calculated, and the flow conductivity and thermal conductivity are corrected respectively with the correction coefficients to obtain the corrected flow conductivity and thermal conductivity of the fourth type of connection pair.

6. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 5, characterized in that, The conductivity of fourth-type connection pairs is corrected using an improved Vermeulen function. The corrected conductivity of fourth-type connection pairs is: ； In the formula, TI ij,final γi represents the corrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid; γ1 and γ2 are empirical coefficients used to fit unsteady mass or heat transfer processes; δi i δi represents the pressure or temperature of the i-th matrix mesh; δ0 represents the initial pressure or temperature of the mesh at the start of the simulation; δi ... j The pressure or temperature of the j-th natural fracture grid; TI ij is the uncorrected flow or thermal conductivity between the i-th matrix grid and the j-th natural fracture grid.

7. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 1, characterized in that, The fully implicit residual equation set includes phase quality residual equations and energy residual equations, specifically: ； Where, when χ=β, it represents the β-phase quality residual equation: ； Where, when χ=T, it represents the energy residual equation: ； In the formula, R χi ρ represents the mass or energy residual in the i-th grid. βi,j λ is the density of phase β at the interface between the i-th and j-th grids; βi,j TI represents the mobility of phase β at the interface between the i-th and j-th grids. ij p represents the flow or conductivity between connected grids. W W represents the bottom hole pressure. i V is the well index of the i-th grid; i Let be the volume of the i-th grid; Δt is the time step.

8. The numerical simulation method for deep fractured reservoir thermal-fluid coupling using a hybrid fracture model according to claim 1, characterized in that, The multiphase flow mass conservation equation includes the flow terms, accumulation terms, and well source-sink terms for each phase of fluid in the corresponding grid system; The energy conservation equation includes convective heat transfer terms, heat conduction terms, heat storage terms, and well-source heat sink terms.

9. A numerical simulation system for coupled thermal flux in deep fractured reservoirs using a hybrid fracture model, characterized in that, include: The basic parameter acquisition module is used to acquire basic parameters of the study area; The mesh construction module is used to establish a computational domain based on the basic parameters, and to construct a matrix mesh system, a natural fracture mesh system, and a large-scale fracture mesh system within the same computational domain; wherein, the large-scale fracture mesh system is explicitly represented using an embedded discrete fracture model, and the natural fracture mesh system and the matrix mesh system are represented using a dual-pore dual-permeability model; The model building module is used to establish the multiphase flow mass conservation equation and energy conservation equation in the matrix grid system, the natural fracture grid system, and the large-scale fracture grid system, respectively. The connectivity calculation module is used to establish four types of mesh connectivity pairs, calculate the flow conductivity and thermal conductivity for each type of mesh connectivity pair, and dynamically correct the conductivity of the fourth type of connectivity pair based on the improved Vermeulen function to obtain the corrected flow conductivity and thermal conductivity. The residual equation system establishment module is used to perform finite volume discretization on the multiphase flow mass conservation equation and energy conservation equation based on the corrected flow conductivity and thermal conductivity, and establish a fully implicit residual equation system. The iterative solution module is used to solve the fully implicit residual equations using the Newton-Raphson iterative method and output the thermal-fluid coupling numerical simulation results of each grid system at different time steps.

Images