Anisotropic medium three-dimensional transient geothermal field forward method, device and medium

Abstract

The application discloses a kind of anisotropic medium three-dimensional transient geothermal field forward method, equipment and medium, by modeling target area, then target area is spatially discretized and time discretization, in space, using the iterative calculation of spatial domain total field temperature of space-wave number hybrid domain iterative method;In time, in combination with the recurrence formula of geothermal field of explicit difference format, the initial spatial domain total field temperature of next time node is calculated, then the spatial domain total field temperature of next time node is calculated using the iterative calculation of spatial domain total field temperature of space-wave number hybrid domain iterative method;In this way, the spatial domain total field temperature of all time nodes is forward calculated.This application realizes the anisotropic medium three-dimensional transient geothermal field forward, can reflect the dynamic change of geothermal field, combine space-wave number hybrid domain iterative method with geothermal field recurrence technique, ensure fine geothermal field forward, reduce the cost of calculation and storage, improve the calculation efficiency.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geothermal field numerical simulation, and particularly relates to a three-dimensional transient geothermal field forward method of anisotropic medium, equipment and medium. BACKGROUND

[0002] In geothermal exploration, due to the universal distribution of non-uniformity of underground medium, studying the anisotropy of the medium is of great significance for the processing and interpretation of geothermal data. The complexity of geological structure and the diversity of thermal reservoir occurrence form lead to the problems of long calculation time and high storage requirement of the existing three-dimensional geothermal field forward method in simulating large-scale complex geological conditions.

[0003] At present, the geothermal field forward method is mainly a steady-state forward method. For example, the patent CN202111423009X discloses a three-dimensional steady-state heat conduction geothermal field forward method, device, equipment and medium. However, the geothermal field is actually dynamic, i.e. transient. Therefore, it is of great practical significance to realize the three-dimensional transient geothermal field forward calculation of anisotropic medium. Although there are some researches on transient geothermal field, the complexity of geological structure and the diversity of thermal reservoir occurrence form lead to the problems of long calculation time and high storage requirement of the existing three-dimensional geothermal field forward method in simulating large-scale complex geological conditions. In view of this problem, the present application provides a fine and efficient three-dimensional transient geothermal field forward technology of anisotropic medium, which provides important technical support for fine inversion of large-scale geothermal data. SUMMARY

[0004] In view of the deficiencies in the prior art, the present application aims to provide a three-dimensional transient geothermal field forward method of anisotropic medium, which realizes fine and efficient three-dimensional transient geothermal field forward of anisotropic medium.

[0005] In the first aspect, a three-dimensional transient geothermal field forward method of anisotropic medium is provided, which comprises:

[0006] S1: selecting a target region containing an anomaly body and constructing a target model;

[0007] S2: performing spatial discretization and time discretization on the target model to obtain a series of sampling nodes and a series of time nodes;

[0008] S3: giving the values of thermal physical parameters on each sampling node, wherein the thermal physical parameters include thermal conductivity, heat generation rate and heat capacity, and the thermal conductivity of the sampling node in the anomaly body region has anisotropy;

[0009] S4: loading initial conditions and boundary conditions;

[0010] S5: Calculate the total temperature of the spatial domain at the current time point using the spatial-wavenumber hybrid domain iterative method;

[0011] S6: Combine the recursive formula of the geothermal field in the explicit difference scheme to calculate the initial spatial domain total field temperature at the next time node, and then use the space-wavenumber hybrid domain iterative method to calculate the spatial domain total field temperature at the next time node; and so on, to perform forward simulation of the spatial domain total field temperature at all time nodes.

[0012] According to the first aspect, in one possible implementation, in step S2, when performing spatial discretization, the target region is spatially divided into a three-dimensional mesh, with uniform division in the x and y directions and uniform or non-uniform division in the z direction.

[0013] According to the first aspect, in one possible implementation, the thermal conductivity at each sampling node includes background thermal conductivity and abnormal thermal conductivity, the heat generation rate at each sampling node includes background heat generation rate and abnormal heat generation rate, and the heat capacity at each sampling node includes background heat capacity and abnormal heat capacity.

[0014] Among them, thermal conductivity For a tensor, the expression is as follows:

[0015]

[0016] In the formula, These are the nine components of thermal conductivity; anomalous thermal conductivity exhibits anisotropy.

[0017] According to the first aspect, in one possible implementation, step S5 includes:

[0018] S51: The temperature of the spatial domain background field is obtained based on the background field control equation and the thermophysical parameters of the background field, and is used as the initial total temperature of the spatial domain.

[0019] S52: Perform a two-dimensional Fourier transform in the horizontal direction on the anomaly field control equations to obtain the space-wavenumber domain anomaly field control equations.

[0020] S53: The temperature of the spatial-wavenumber domain anomaly field is calculated based on the total temperature of the spatial domain and the thermophysical parameter values ​​of the anomaly field in the current round, combined with the control equation of the spatial-wavenumber domain anomaly field.

[0021] S54: Perform a two-dimensional inverse Fourier transform on the temperature of the anomalous field in the space-wavenumber domain to obtain the temperature of the anomalous field in the space domain;

[0022] S55: The new total field temperature of the spatial domain is obtained based on the background field temperature and the anomalous field temperature of the spatial domain.

[0023] S56: judging whether the new spatial domain total field temperature meets an iterative convergence condition, if yes, outputting the new spatial domain total field temperature as the spatial domain total field temperature of the current time node; otherwise, taking the new spatial domain total field temperature as the current round spatial domain total field temperature in the next round of iteration, and returning to step S53.

[0024] According to the first aspect, in a possible implementation manner, the background field control equation is represented as follows:

[0025]

[0026] wherein, is a background field temperature, is a background thermal conductivity, is a background heat generation rate, is a background heat capacity, and t is time;

[0027] The spatial-wave number domain anomaly field control equation is represented as follows:

[0028]

[0029] wherein, is an anomaly thermal conductivity, is an anomaly heat capacity, is a spatial-wave number domain anomaly field temperature, is a spatial-wave number domain anomaly heat generation rate, are respectively wave numbers in x, y and z directions, is a two-dimensional Fourier transform symbol, i is an imaginary unit, and T is a total field temperature, is a fluid heat capacity, and v x , v y , v z are respectively flow velocities of the fluid in x, y and z directions.

[0030] According to the first aspect, in a possible implementation manner, the iterative convergence condition is that an error is less than a preset value; and an expression of the error is as follows:

[0031]

[0032] wherein, e is the error, are respectively numbers of sampling nodes in x, y and z directions, , are respectively the spatial domain total field temperatures in the n th iteration and the n+1 th iteration.

[0033] According to the first aspect, in a possible implementation manner, the step S6 comprises:

[0034] S61: Combining the recursive formula of the geothermal field in the explicit difference scheme, the spatial-wavenumber domain anomaly temperature of the next time node is obtained by solving the spatial-wavenumber domain anomaly temperature in the last iteration of the current time node using the chasing method.

[0035] S62: Perform a two-dimensional inverse Fourier transform on the spatial-wavenumber domain anomaly temperature at the next time node to obtain the spatial domain anomaly temperature at the next time node.

[0036] S63: The initial total field temperature of the spatial domain at the next time node is obtained based on the spatial domain anomaly field temperature and the spatial domain background field temperature at the next time node;

[0037] S64: Calculate the total spatial field temperature at the next time node using the spatial-wavenumber hybrid domain iterative method;

[0038] S65: Repeat steps S61 to S64 until the total field temperature of the spatial domain at all time points is obtained through forward modeling.

[0039] According to the first aspect, in one possible implementation, the recursive formula for the geothermal field is expressed as follows:

[0040]

[0041] In the formula, The previous time point Space-wavenumber domain anomalous field temperature, For the next time node Space-wavenumber domain anomalous field temperature, Let K be the time interval between adjacent time nodes, G be the coefficient matrix related to the spatial-wavenumber domain anomaly field, G be the coefficient matrix related to the time derivative of the spatial-wavenumber domain anomaly field, and P be the source term.

[0042] According to the first aspect, in one possible implementation, the total field temperature of the spatial domain between adjacent time points is obtained by the following method:

[0043] Get adjacent time nodes and time nodes Spatial domain background field temperature and and obtaining time nodes and time nodes Spatial-wavenumber domain anomalous field temperature and ;

[0044] Calculate adjacent time nodes using interpolation algorithms and time nodes Spatial background field temperature at time t′ and spatial-wavenumber domain anomaly field temperature ;

[0045]

[0046] wherein a is a coefficient constant, ;

[0047] performing two-dimensional Fourier inverse transform on the spatial-wavenumber domain anomaly field temperature at the time t' to obtain a spatial domain anomaly field temperature at the time t';

[0048] adding the spatial domain background field temperature at the time t' and the spatial domain anomaly field temperature to obtain a spatial domain total field temperature at the time t'.

[0049] According to the first aspect, in a possible implementation manner, time intervals between adjacent time nodes in the series of time nodes are equal, or are not equal, or are partially equal.

[0050] The second aspect provides an electronic device, comprising:

[0051] a memory storing a computer program;

[0052] a processor configured to load and execute the computer program to implement steps of the anisotropic medium three-dimensional transient geothermal field forward method.

[0053] The third aspect provides a computer readable storage medium storing a computer program, which, when executed by a processor, implements steps of the anisotropic medium three-dimensional transient geothermal field forward method.

[0054] The present application provides an anisotropic medium three-dimensional transient geothermal field forward method, device and medium, which discretizes a target region in space and time, iteratively calculates a spatial domain total field temperature in space by using a spatial-wavenumber hybrid domain iterative method, combines a geothermal field recursive formula of an explicit difference format to calculate an initial spatial domain total field temperature of a next time node, and then calculates the spatial domain total field temperature of the next time node by using the spatial-wavenumber hybrid domain iterative method, and iteratively calculates spatial domain total field temperatures of all time nodes. BRIEF DESCRIPTION OF DRAWINGS

[0055] In order to more clearly illustrate the technical solutions of the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings described below are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained from these drawings without creative labor.

[0056] Figure 1 is a flow chart of the three-dimensional geothermal field dynamic numerical simulation method provided by the embodiments of the present application;

[0057] Figure 2 is a flow chart of the spatial domain total field temperature at the current time node calculated by the spatial-wave number hybrid domain iterative method provided by the embodiments of the present application;

[0058] Figure 3 is a flow chart of the forward of the spatial domain total field temperature of all time nodes provided by the embodiments of the present application;

[0059] Figure 4 is a target model profile projection diagram provided by the embodiments of the present application, wherein (a) is a target model XOY profile projection diagram, and (b) is a target model XOZ profile projection diagram;

[0060] Figure 5 is a schematic diagram of the calculation results of the present application method and COMSOL Multiphysics software at a sampling node provided by the embodiments of the present application;

[0061] Figure 6 is a schematic diagram of the relative error of the calculation results of the present application method and COMSOL Multiphysics software at a sampling node provided by the embodiments of the present application;

[0062] Figure 7 is a schematic diagram of the calculation results of the present application method and COMSOL Multiphysics software at another sampling node provided by the embodiments of the present application;

[0063] Figure 8 is a schematic diagram of the relative error of the calculation results of the present application method and COMSOL Multiphysics software at another sampling node provided by the embodiments of the present application. DETAILED DESCRIPTION

[0064] In order to make the purpose, technical solutions and advantages of the present application more clear, the technical solutions of the present application will be described in detail below. Obviously, the described embodiments are only some embodiments of the present application, not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of the present application.

[0065] This invention proposes a three-dimensional dynamic numerical simulation method for geothermal fields, such as... Figure 1 As shown, it includes the following steps:

[0066] S1: Select the target region containing the anomaly and construct the target model.

[0067] S2: Spatial and temporal discretization of the target model is performed to obtain a series of sampling nodes and a series of time nodes.

[0068] In this embodiment, during spatial discretization, the target region is spatially divided into a three-dimensional mesh. The x and y directions are uniformly divided, and the z direction is either uniformly or non-uniformly divided; the number of sampling nodes in the x direction is... The number of sampling nodes in the y-direction is The number of sampling nodes in the z-direction is When performing time discretization, a series of time nodes are selected, with the number of time nodes being... .

[0069] S3: Given the thermophysical parameter values ​​at each sampling node, the thermophysical parameters include thermal conductivity, heat generation rate and heat capacity, wherein the thermal conductivity of the sampling nodes in the anomalous body region is anisotropic.

[0070] Specifically, the thermal conductivity at each sampling node includes background thermal conductivity and abnormal thermal conductivity, the heat generation rate at each sampling node includes background heat generation rate and abnormal heat generation rate, and the heat capacity at each sampling node includes background heat capacity and abnormal heat capacity; for sampling nodes without anomalies, the abnormal thermal conductivity, abnormal heat generation rate, and abnormal heat capacity are all 0; the background thermal conductivity, heat generation rate, and heat capacity at sampling nodes with anomalies are set as parameters of a uniform layered model.

[0071] Among them, thermal conductivity For a tensor, the expression is as follows:

[0072]

[0073] In the formula, These are the nine components of thermal conductivity; anomalous thermal conductivity exhibits anisotropy.

[0074] S4: Load initial and boundary conditions.

[0075] Specifically, boundary conditions refer to known temperature, heat flux density, or heat transfer coefficient at the boundary. Initial conditions are given at an initial time. The ground temperature field value at that time, the initial temperature of each sampling node in the target area can be given as raw data, or it can be calculated based on boundary conditions and model parameters.

[0076] S5: Calculate the spatial-domain total-field temperature at the current time node by using the spatial-wave-number hybrid domain iterative method. Specifically, as shown in FIG. 5, the step S5 comprises: Figure 2

[0077] S51: Obtain the spatial-domain background-field temperature based on the background-field governing equation and the values of the thermophysical parameters of the background field (i.e. the values of the background thermal conductivity, the background heat generation rate and the background heat capacity), and take it as the initial spatial-domain total-field temperature.

[0078] The background-field governing equation is expressed as follows:

[0079]

[0080] wherein, is the background-field temperature, is the background thermal conductivity, is the background heat generation rate, is the background heat capacity, and t is time.

[0081] S52: Perform a two-dimensional Fourier transform in the horizontal direction on the anomalous-field governing equation to obtain the spatial-wave-number domain anomalous-field governing equation.

[0082] Specifically, the anomalous-field governing equation is expressed as follows:

[0083]

[0084] wherein, is the anomalous-field temperature, is the background thermal conductivity, is the background heat capacity, is the anomalous thermal conductivity, is the anomalous heat generation rate, is the anomalous heat capacity, and T is the total-field temperature, is the fluid heat capacity, and v is the flow rate of the fluid.

[0085] Perform a two-dimensional Fourier transform in the horizontal direction on the anomalous-field governing equation to obtain the spatial-wave-number domain anomalous-field governing equation (i.e. the one-dimensional spatial-domain anomalous-field governing equation under two-dimensional wave-number domain), which is expressed as follows:

[0086]

[0087] wherein, is the spatial-wave-number domain anomalous-field temperature, is the spatial-wave-number domain anomalous heat generation rate, are respectively the wave numbers in the x and y directions, is the two-dimensional Fourier transform symbol, i is the imaginary unit, and v x , v​y , v z are the fluid velocities in the x, y, z directions, respectively.

[0088] S53: According to the current round space domain total field temperature, the abnormal field thermophysical parameter value (i.e. the values of abnormal thermal conductivity, abnormal heat generation rate and abnormal heat capacity), combined with the space-wave number domain abnormal field control equation, the one-dimensional finite element method based on quadratic interpolation is used to solve and calculate the space-wave number domain abnormal field temperature.

[0089] S54: The two-dimensional Fourier inverse transform is performed on the space-wave number domain abnormal field temperature to obtain the space domain abnormal field temperature.

[0090] S55: Based on the space domain background field temperature and the space domain abnormal field temperature, a new space domain total field temperature is obtained.

[0091] S56: It is judged whether the new space domain total field temperature meets the iterative convergence condition or not. If yes, the new space domain total field temperature is output as the space domain total field temperature of the current time node; otherwise, the new space domain total field temperature is taken as the current round space domain total field temperature in the next iteration, and the step S53 is returned.

[0092] In this embodiment, the iterative convergence condition is that the error is less than a preset value; and the expression of the error is as follows:

[0093]

[0094] In the formula, e is the error, are the number of sampling nodes in the x, y, z directions, respectively, , are the space domain total field temperatures of the n th and n+1 th iterations, respectively. That is, if the error e between the space domain total field temperature obtained in the current iteration round and the space domain total field temperature obtained in the last iteration round is less than a set value, the convergence condition is met, the iteration is ended, and the space domain total field temperature is output; otherwise, the convergence condition is not met, and the steps S6-S9 are repeated to continue the iteration calculation until the convergence is achieved.

[0095] S6: The initial space domain total field temperature of the next time node is calculated by combining the explicit difference format of the geothermal field recursive formula, and then the space-wave number hybrid domain iteration method is used to calculate the space domain total field temperature of the next time node. In this way, the space domain total field temperatures of all time nodes are calculated.

[0096] Specifically, as shown in the figure, the step S6 includes: Figure 3

[0097] ​S61: Combining the recursive formula of the geothermal field in the explicit difference scheme, the spatial-wavenumber domain anomaly temperature of the next time node is obtained by solving the spatial-wavenumber domain anomaly temperature in the last iteration of the current time node using the chasing method.

[0098] The recursive formula for the geothermal field is expressed as follows:

[0099]

[0100] In the formula, The previous time point Space-wavenumber domain anomalous field temperature, For the next time node The temperature of the spatial-wavenumber domain anomaly field, where K is the coefficient matrix related to the spatial-wavenumber domain anomaly field. Let K be the time interval between adjacent time points, G be the coefficient matrix related to the time derivative of the space-wavenumber domain anomaly field, and P be the source term. In the process of solving the control equations of the space-wavenumber domain anomaly field using the finite element method to obtain the temperature of the space-wavenumber domain anomaly field, the values ​​of K, G, and P can be obtained through element analysis and overall synthesis. The above equation can be simplified to a diagonal system of equations. In this system of diagonal equations, Q is a pentagonal matrix, X is the unknown variable to be solved, and B is the right-hand side term. The chasing method is used to solve this system of diagonal equations quickly. The chasing method for solving ordinary differential equations has high parallelism, which greatly reduces computation time and memory usage.

[0101] S62: Perform a two-dimensional inverse Fourier transform on the spatial-wavenumber domain anomaly temperature at the next time node to obtain the spatial domain anomaly temperature at the next time node.

[0102] S63: Add the spatial domain anomaly field temperature to the spatial domain background field temperature at the next time node to obtain the initial spatial domain total field temperature at the next time node.

[0103] S64: Calculate the total spatial field temperature at the next time node using the spatial-wavenumber hybrid domain iterative method (corresponding to the aforementioned steps S53~S56);

[0104] S65: Repeat steps S61 to S64 until the total field temperature of the spatial domain at all time points is obtained through forward modeling.

[0105] In a preferred embodiment of the present invention, it further includes:

[0106] S7: Use an interpolation algorithm to obtain the total field temperature of the spatial domain between adjacent time nodes.

[0107] Step S7 specifically includes:

[0108] S71: Obtain adjacent time nodes and time nodes spatial domain background field temperature and , and the acquisition time nodes and time nodes spatial-wave number domain anomaly field temperature and ;

[0109] S72: calculating the spatial domain background field temperature at the moment t' between adjacent time nodes and time nodes and the spatial-wave number domain anomaly field temperature and spatial-wave number domain anomaly field temperature ;

[0110]

[0111] wherein a is a coefficient constant, ;

[0112] S73: performing two-dimensional Fourier inverse transform on the spatial-wave number domain anomaly field temperature at the moment t' to obtain the spatial domain anomaly field temperature at the moment t';

[0113] S74: adding the spatial domain background field temperature at the moment t' and the spatial domain anomaly field temperature to obtain the spatial domain total field temperature at the moment t'.

[0114] It should be noted that the time intervals between adjacent time nodes in the series of time nodes can be equal, unequal, or partially equal.

[0115] The embodiment of the present application further provides an electronic device, comprising:

[0116] a memory storing a computer program;

[0117] a processor configured to load and execute the computer program to implement the steps of the anisotropic medium three-dimensional transient geothermal field forward method as described in the above embodiment.

[0118] The embodiment of the present application further provides a computer readable storage medium storing a computer program, wherein the computer program is executed by a processor to implement the steps of the anisotropic medium three-dimensional transient geothermal field forward method as described in the above embodiment.

[0119] The effect of the scheme provided by the present application is verified below in combination with a specific embodiment.

[0120] The computer configuration for testing is Intel(R) Core(TM) i7-11800H, main frequency 2.30 GHz, memory 16 GB, 64-bit win11 system.

[0121] The target model XOY plane and XOZ profile projection are as shown in Figure 4 Figs. (a) and (b), the calculation range is x direction -5000~5000 m, y direction -5000~5000 m, and z direction 0~10000 m.

[0122] The background is a full-space medium, the background thermal conductivity is 2 W / (m•℃), the heat capacity is 5×10 5 J / (m 3 •℃), the anisotropic abnormal thermal conductivity components are 3, 1, 3, 1, 3, 2, 3, 2, 4 W / (m•℃), the abnormal heat capacity is 1×10 6 J / (m3•℃), the background has no endogenous heat, and the abnormal heat generation rate is 4×10 -6 W / m 3 . The upper boundary adopts a first type of boundary condition (temperature value is known), the boundary initial temperature value is 10℃, and the lower boundary adopts a second type of boundary condition (heat flow density value is known), the heat flow density value is 41.86 mW / m 2 . The spatial grid is divided into 51×51×51, the time node is 21, and the iteration convergence precision (preset value) is 10 -4 . The forward is performed by using the method and the COMSOL Multiphysics software respectively, and the temperature change results with time at sampling nodes and sampling nodes are selected. Figure 5 The calculation results of the two methods at sampling nodes , Figure 6 the relative error diagram, Figure 7 the calculation results of the two methods at sampling nodes , Figure 8 and the relative error diagram. It can be seen that the relative errors of the calculation results of the two methods are less than 0.06%. Under the same division, the time consumption of the method is 58 s, the memory occupation is 0.17 GB, the calculation time of the COMSOL Multiphysics is 287 s, and the memory occupation is 7.79 GB, which shows the advantages of the method.

[0123] It can be understood that the same or similar parts in the above embodiments can be mutually referred to, and the contents not described in detail in some embodiments can be referred to the same or similar contents in other embodiments.

[0124] Those skilled in the art will appreciate that embodiments of the present application can be readily used as a method, a system or a computer program product. Accordingly, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present application can take the form of a computer program product on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROMs, optical storage devices, etc.) embodying computer readable program code.

[0125] The present application is described in reference to the flowchart and / or block diagrams of the methods, apparatus (systems) and computer program products according to embodiments of the present application. It will be understood that each block of the flowchart and / or block diagrams, and combinations of blocks in the flowchart and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general purpose computer, special purpose computer, embedded processing system or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks.

[0126] These computer program instructions can also be stored in a computer- readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instructions which implement the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks.

[0127] These computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks. Figure 1 one or more functions specified in the flowchart and / or block diagram block or blocks.

[0128] Although the embodiments of the present application have been described above, it is understood that the above embodiments are merely exemplary and are not to be construed as limiting the present application, and that modifications, substitutions, replacements and variations of the above embodiments can be made by those skilled in the art within the scope of the present application.

Claims

1. A forward modeling method for three-dimensional transient geothermal field in anisotropic media, characterized in that, include: S1: Select the target region containing the anomaly and construct the target model; S2: Spatial and temporal discretization of the target model is performed to obtain a series of sampling nodes and a series of time nodes; S3: Given the thermophysical parameter values ​​at each sampling node, the thermophysical parameters include thermal conductivity, heat generation rate and heat capacity, wherein the thermal conductivity of the sampling nodes in the anomalous body region is anisotropic; S4: Load initial and boundary conditions; S5: Calculate the total temperature of the spatial domain at the current time point using the spatial-wavenumber hybrid domain iterative method; S6: Combine the recursive formula of the geothermal field in the explicit difference scheme to calculate the initial spatial domain total field temperature at the next time node, and then use the spatial-wavenumber hybrid domain iterative method to calculate the spatial domain total field temperature at the next time node. And so on, the total field temperature of the spatial domain at all time points is calculated. Step S6 includes: S61: Combining the recursive formula of the geothermal field in the explicit difference scheme, the spatial-wavenumber domain anomaly temperature of the next time node is obtained by solving the spatial-wavenumber domain anomaly temperature in the last iteration of the current time node using the chasing method. S62: Perform a two-dimensional inverse Fourier transform on the spatial-wavenumber domain anomaly temperature at the next time node to obtain the spatial domain anomaly temperature at the next time node. S63: The initial total field temperature of the spatial domain at the next time node is obtained based on the spatial domain anomaly field temperature and the spatial domain background field temperature at the next time node; S64: Calculate the total spatial field temperature at the next time node using the spatial-wavenumber hybrid domain iterative method; S65: Repeat steps S61 to S64 until the total field temperature of the spatial domain at all time points is obtained through forward modeling.

2. The method for forward modeling the three-dimensional transient geothermal field of anisotropic media according to claim 1, characterized in that, The thermal conductivity at each sampling node includes background thermal conductivity and abnormal thermal conductivity; the heat generation rate at each sampling node includes background heat generation rate and abnormal heat generation rate; and the heat capacity at each sampling node includes background heat capacity and abnormal heat capacity. Among them, thermal conductivity For a tensor, the expression is as follows: ； In the formula, These are the nine components of thermal conductivity; anomalous thermal conductivity exhibits anisotropy.

3. The forward modeling method for three-dimensional transient geothermal field of anisotropic media according to claim 1, characterized in that, Step S5 includes: S51: The temperature of the spatial domain background field is obtained based on the background field control equation and the thermophysical parameters of the background field, and is used as the initial total temperature of the spatial domain. S52: Perform a two-dimensional Fourier transform in the horizontal direction on the anomaly field control equations to obtain the space-wavenumber domain anomaly field control equations. S53: The temperature of the spatial-wavenumber domain anomaly field is calculated based on the total temperature of the spatial domain and the thermophysical parameter values ​​of the anomaly field in the current round, combined with the control equation of the spatial-wavenumber domain anomaly field. S54: Perform a two-dimensional inverse Fourier transform on the temperature of the anomalous field in the space-wavenumber domain to obtain the temperature of the anomalous field in the space domain; S55: The new total field temperature of the spatial domain is obtained based on the background field temperature and the anomalous field temperature of the spatial domain. S56: Determine whether the new spatial domain total field temperature satisfies the iteration convergence condition. If it does, output the new spatial domain total field temperature as the spatial domain total field temperature at the current time node; otherwise, use the new spatial domain total field temperature as the spatial domain total field temperature of the current round in the next iteration and return to step S53.

4. The forward modeling method for three-dimensional transient geothermal field of anisotropic media according to claim 3, characterized in that, The background field governing equations are expressed as follows: ； In the formula, For the background field temperature, Background thermal conductivity, For background heat generation rate, The background heat capacity is t, and the time is t. The governing equations for the space-wavenumber domain anomaly field are expressed as follows: ； In the formula, For abnormal thermal conductivity, For abnormal heat capacity, For the temperature of the anomalous field in the space-wavenumber domain, For the space-wavenumber domain anomalous heat generation rate, They are respectively Wave number in direction, This is the symbol for a two-dimensional Fourier transform, where i is the imaginary unit and T is the total field temperature. For the fluid heat capacity, v x v y v z These represent the fluid velocities in the x, y, and z directions, respectively.

5. The forward modeling method for three-dimensional transient geothermal field of anisotropic media according to claim 3, characterized in that, The iteration convergence condition is that the error is less than a preset value; the expression for the error is as follows: ； In the formula, e represents the error. These represent the number of sampling nodes in the x, y, and z directions, respectively. , These are the total field temperatures of the spatial domain in the nth and (n+1)th iterations, respectively.

6. The method for forward modeling the three-dimensional transient geothermal field of anisotropic media according to claim 1, characterized in that, The recursive formula for the geothermal field is expressed as follows: ； In the formula, The previous time point Space-wavenumber domain anomalous field temperature, For the next time node Space-wavenumber domain anomalous field temperature, Let K be the time interval between adjacent time nodes, G be the coefficient matrix related to the spatial-wavenumber domain anomaly field, G be the coefficient matrix related to the time derivative of the spatial-wavenumber domain anomaly field, and P be the source term.

7. The method for forward modeling the three-dimensional transient geothermal field of anisotropic media according to claim 3, 4, or 5, characterized in that, The total field temperature in the spatial domain between adjacent time points is obtained using the following method: Get adjacent time nodes and time nodes Spatial domain background field temperature and and obtaining time nodes and time nodes Spatial-wavenumber domain anomalous field temperature and ; Calculate adjacent time nodes using interpolation algorithms and time nodes Spatial background field temperature at time t′ Temperature of the space-wavenumber domain anomaly field ; ； In the formula, 'a' is a constant coefficient. ; A two-dimensional inverse Fourier transform is performed on the spatial-wavenumber domain anomaly temperature at time t′ to obtain the spatial domain anomaly temperature at time t′. The total temperature of the spatial domain at time t′ is obtained by adding the background temperature of the spatial domain to the anomalous temperature of the spatial domain.

8. An electronic device, characterized in that, include: A memory that stores computer programs; A processor for loading and executing the computer program to implement the steps of the forward modeling method for three-dimensional transient geothermal fields of anisotropic media as described in any one of claims 1 to 7.

9. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the steps of the forward modeling method for the three-dimensional transient geothermal field of anisotropic media as described in any one of claims 1 to 7.

Images