Three-dimensional transient prediction method, device and electronic equipment for thermal conductivity of deep formations

Abstract

The present invention provides a three-dimensional transient prediction method, device and electronic device for thermal conductivity of deep formations. The method uses underground heat flow values, thermal physical parameters of each formation and fluid migration law as prior information, conducts transient finite element temperature numerical simulation and constructs a regularized objective function, constructs a Hessian matrix and a gradient vector based on the Gauss-Newton algorithm, uses the Jacobian-free Krylov subspace technology based on the time chain derivation law to solve the product of the Jacobian matrix and an arbitrary vector to avoid solving and storing large dense Jacobian matrices and obtain model corrections, searches for the optimal model step size and updates model parameters based on the Wolfe criterion, and cyclically predicts so that the data fitting error is less than the preset fitting error, and outputs the current prediction model parameters as the optimal thermal conductivity prediction result. The present invention can intuitively characterize the distribution characteristics of the thermal conductivity of deep formation media in the time domain, has high prediction accuracy and strong practicality.

Description

Technical Field

[0001] The present invention relates to the field of thermal conductivity prediction, and in particular to a three-dimensional transient prediction method, device and electronic equipment for thermal conductivity of deep formations. Background Art

[0002] As a new type of energy, geothermal resources are attracting more and more attention due to their unique advantages that other energy sources cannot match, such as low cost, sustainable use and environmental protection. Vigorously promoting the development and utilization of geothermal resources and improving the energy structure are of great significance to solving the increasingly serious global environmental problems and achieving the "dual carbon" goals. Thermal conductivity, thermal diffusion coefficient, specific heat capacity, etc. are very important rock thermophysical parameters in deep strata. They are the basis for quantitative evaluation of temperature changes. Accurate calculation of temperature is very important for clarifying underground heat transfer patterns, geothermal genesis mechanisms, and delineating geothermal resource target areas.

[0003] Thermal conductivity is one of the most important parameters in the thermophysical properties of rocks. Due to the different geological structures and evolutionary processes inside the earth, the rock properties are also very different, resulting in obvious heterogeneity in the spatial distribution of the thermal conductivity of rocks, which in turn causes changes in the geothermal characteristics of the region. At present, the prediction of rock thermal conductivity is usually achieved through laboratory measurements of drill cuttings or core samples or through geophysical logging. However, these methods only characterize the steady-state thermal conductivity of rocks at point and line scales, ignoring the impact of the flow of underground thermal fluids on the thermal conductivity of rocks, and limiting the accurate prediction of regional time-domain rock thermal conductivity. This has strong limitations for research in the fields of high-temperature geothermal target site selection and thermal genesis analysis in the study area. Summary of the invention

[0004] In order to solve the technical problem that the existing rock thermal conductivity prediction method can only characterize the steady-state rock thermal conductivity at the point and line scales and ignores the influence of the flow of underground thermal fluid on the rock thermal conductivity, the present invention proposes a three-dimensional transient prediction method, device and electronic device for deep formation thermal conductivity. The method uses underground heat flow values, thermal physical parameters of each formation and fluid migration law as prior information, carries out transient finite element temperature numerical simulation and constructs a regularized objective function, constructs the Hessian matrix and gradient vector based on the Gauss-Newton algorithm, uses the Jacobian-free Krylov subspace technology based on the time chain derivation law to solve the product of the Jacobian matrix and any vector to avoid solving and storing large dense Jacobian matrices and obtain model corrections, searches for the optimal model step size to update model parameters using the Wolfe criterion, and cyclically predicts so that the data fitting error is less than the preset fitting error, outputs the current prediction model parameters as the optimal thermal conductivity prediction result, and realizes the three-dimensional transient prediction of deep medium thermal conductivity.

[0005] According to a first aspect of the present invention, a method for three-dimensional transient prediction of thermal conductivity of deep formations is provided, comprising the following steps:

[0006] Obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly in the study area and give boundary conditions, conduct transient finite element temperature numerical simulation based on the distribution of thermophysical parameters in the study area and the boundary conditions, and obtain the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively;

[0007] According to the underground space temperature field d obs , using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm;

[0008] According to the prediction model correction amount Δm, an optimal prediction model search step length satisfying the Wolfe criterion is obtained, and prediction model parameters are updated;

[0009] According to the updated prediction model parameters, transient finite element temperature numerical simulation is carried out to obtain the underground space temperature field under the model parameters, and the data fitting difference between the measured data is calculated. If it is less than the preset fitting difference condition, the current prediction model parameters are the optimal thermal conductivity prediction results; otherwise, return to the prediction method and solve the new prediction model correction value Δm.

[0010] Preferably, the groundwater migration law and its response mechanism to thermal conductivity include:

[0011] Assume that the groundwater movement rates in the x, y, and z directions of the study area are v x 、v y 、v z , the fluid migrates and immerses into the surrounding medium, causing the thermal conductivity of the medium to change, satisfying:

[0012] k=(1-φ)·k s +φ·k f

[0013] Among them, k, k s , k fare the true thermal conductivity, medium thermal conductivity and fluid thermal conductivity, respectively, with the unit of W / m·K; φ is the medium porosity, dimensionless.

[0014] Preferably, the boundary conditions satisfy:

[0015]

[0016] Among them, α and β are constants, α=1, β=T up When is the upper boundary condition; α=0, β=q down When α=0, β=0, it is the lower boundary condition; when α=0, β=0, it is the surrounding boundary condition. T is the temperature, unit ℃, T up is the surface temperature, unit is °C, n is the direction of the boundary normal, dimensionless, k is the thermal conductivity, unit is W / m·K, q down is the heat flow value, unit is W / m 2 .

[0017] Preferably, the transient finite element temperature numerical simulation adopts a hexahedral structure discretization method, and the temperature T, thermal conductivity k, and fluid heat capacity μ in the hexahedral unit are f 、Medium heat capacity μ s , fluid migration rate v, and heat source Q are all linearly interpolated:

[0018]

[0019]

[0020] Where i is the node number, T i , k i , μ fi , μ si 、v i , Q i They are the temperature, thermal conductivity, fluid heat capacity, medium heat capacity, fluid flow rate, and heat source of each node in the unit, respectively, in units of ℃, W / m·K, J / (m 3 ·K)、J / (m 3 ·K), m / s, W / m 3 ; N i is a shape function that satisfies:

[0021]

[0022] Among them, ξ i , η i , i is the coordinate of node i on the parent unit, dimensionless; the relationship formula between ξ, η, ζ in the parent unit and the coordinates x, y, z in the child unit is:

[0023]

[0024] Among them, x 0 ,y 0 、z 0 are the midpoints of the subunit in the x, y, and z directions respectively; a, b, and c are the side lengths of the subunit in the x, y, and z directions respectively, in m.

[0025] Preferably, the step of carrying out transient finite element temperature numerical simulation in combination with the distribution of thermophysical property parameters in the study area and the boundary conditions to obtain the temperature field of the underground space at a given time includes:

[0026] Given the transient heat transfer underground space temperature field differential equation is:

[0027]

[0028] Where k is the thermal conductivity, unit W / m·K, T is the temperature, unit ℃, μ f is the heat capacity of the fluid, in J / (m 3 ·K), μ s is the heat capacity of the medium, in J / (m 3 K), v is the fluid velocity, unit is m / s, Q is the heat source, unit is W / m 3 , t is time, unit is s, is the divergence, is the gradient.

[0029] The differential equation of the underground space temperature field is multiplied by the impulse function δT and integrated:

[0030]

[0031] Among them, Ω is the study area, Γ is the boundary, and δT is the impulse function.

[0032] According to the Hamiltonian operator operation rules and Ogaard formula in field theory, and substituting the boundary conditions, the first term in the above formula is rewritten as:

[0033]

[0034] Therefore, the integral equation of the temperature field of the underground space with transient heat transfer is:

[0035]

[0036] The integral equation of the underground space temperature field is discretized and linearly interpolated to obtain the unit integral of each term in the integral equation:

[0037]

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] Where t is time, unit is s, δT e is the impulse function array in the unit, T e is the temperature array in the unit, K 1e ,K 2e ,K 3e ,G e is the stiffness matrix within the element, P 1e ,P 2e is the source vector in the element, and the stiffness matrix and the coefficients in the source vector are obtained by solving:

[0044]

[0045]

[0046]

[0047]

[0048]

[0049]

[0050]

[0051] Among them, k 1ij ,k 2ij ,k 3ij ,g ij ,p 1ij ,p 2ij The stiffness matrix K is 1e ,K 2e ,K 3e ,G e , source vector P 1e ,P 2e The elements in the , i, j, l are node labels, a, b, c are the side lengths of the subunit in the x, y, z directions, respectively, in m, ξ i , η i , i is the coordinate of node i on the parent unit.

[0052] The numbers of each node of each unit are extended to the element positions corresponding to the study area, and the linear equations of the synthesized overall stiffness matrix and source vector are obtained:

[0053]

[0054] Among them, δT is the impulse function array of the study area, T is the temperature array in the study area, K 1 ,K 2 ,K 3 , G are the stiffness matrices in the study area, P 1 ,P 2 is the source vector in the study area, ignoring δT T , after simplification:

[0055]

[0056] Where K = K 1 +K 2 +K 3 , P=P 1 +P 2 .

[0057] The temperature field of the underground space in the transient underground heat transfer is differentiated at time t. After Δt, the temperature field is changed from T 1 Change to T 2 ,but In addition, the interpolated temperature within the Δt time Substituting T in the above formula, the transient temperature field linear equations are:

[0058]

[0059] By solving the above equation and realizing finite element numerical simulation, we can get the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively.

[0060] Preferably, the objective function of the transient temperature forward modeling constructed by Tikhonov regularization is:

[0061] Φ=||D[d(m)-d obs ]|| 2 +λ||W(mm apr )|| 2

[0062] Where D is the data weighting matrix; d obs is the observed data, i.e., the underground space temperature field T(x, y, z); m is the current prediction model parameter vector; d(m) is the simulated data under the current prediction model parameters; W is the smooth constraint matrix; λ is the damping factor; mapr Prior information for the prediction model;

[0063] Based on the linear equations of the overall stiffness matrix and source vector in the finite element, the thermal conductivity k in the linear equations is derived as follows:

[0064]

[0065]

[0066]

[0067] Among them, K, G are the stiffness matrices in the study area, P 1 ,P 2 is the source vector in the study area, Δt is the time interval in seconds, k is the thermal conductivity array in W / m·K, T i , T i+1 are the temperature arrays at time i and i+1 respectively. At the initial time (t=0), the above formula is rewritten as:

[0068]

[0069] Among them, T 0 is the temperature array at the initial moment. To ensure that the temperature value and thermal conductivity are non-negative, the data parameters and model parameters are set as the logarithm of temperature and the logarithm of thermal conductivity during the prediction process: d = lnT, m = lnk; according to the chain rule, the elements in the Jacobian matrix are:

[0070]

[0071] Where i = 1, 2, ..., N, N is the number of data parameters, j = 1, 2, ..., M, M is the number of model parameters, d i is the data parameter of the ith node in the study area, m j represents the model parameters of the jth node in the study area, T i represents the temperature of the ith node in the study area, k j represents the thermal conductivity of the jth node in the study area, then the Jacobian matrix is ​​written as:

[0072] J=RGR'

[0073] Among them, R and R' are diagonal matrices, R'=diag(k 1 ,k 2 ,…,k M ), G is a matrix composed of elements G(i,j);

[0074] For any vector consisting of dimensionless elements v = (v 1 ,v 2 ,…,v M ) T , and the conversion formula based on the finite element stiffness matrix is:

[0075]

[0076] Among them, K(v) is the stiffness matrix when the thermal conductivity array takes v; and when t = 0, there is

[0077]

[0078] Among them, e 1 ,e 2 ,...,e M are all unit vectors, so the product of the Jacobian matrix J and any vector v can be found.

[0079] Similarly, the transpose of the Jacobian matrix is ​​calculated and the arbitrary vector y=(y 1 ,y 2 ,…,y N ) T The product of

[0080] According to the Gauss-Newton algorithm, ignoring the second-order information term, the following linear equations are obtained to solve the model correction value Δm:

[0081] HΔm=-g

[0082] Where H is the Hessian matrix, H = (J T D T DJ+λW T W), g is the gradient vector, g = [J T D T D(dd obs )+λW T W(mm apr )].

[0083] Preferably, obtaining the optimal prediction model search step length satisfying the Wolfe criterion according to the prediction model correction amount Δm needs to satisfy the following inequality:

[0084]

[0085] Where α is the search step length, dimensionless, and c is a dimensionless number. In the process of obtaining the search step length, let α be 1 as the initial value. If the inequality is not satisfied, let

[0086] α=ρα

[0087] Among them, ρ is a decimal greater than 0 and less than 1, and the cycle is repeated until the Wolfe criterion is met; therefore, the update formula of the model parameters is:

[0088] m k+1 =m k +αΔm k

[0089] Among them, m k 、m k+1 are the model parameters at the kth and k+1th iterations, Δm k is the model correction amount obtained in the kth iteration.

[0090] Preferably, the calculation formula of the data fitting error is as follows:

[0091]

[0092] RMS is the root mean square, which is used to measure the deviation between the measured value and the predicted value. The closer it is to 0, the smaller the deviation between the measured value and the predicted value, indicating that the prediction result is more reliable; d(m) i d obs i are the i-th predicted value and measured value respectively, n is the number of observed data, dimensionless;

[0093] If the data fitting error is less than the preset fitting error, the iteration is stopped and the current prediction model parameters are output as the optimal prediction result; otherwise, the underground space temperature field d(m) is updated according to the model parameters and the new prediction model correction Δm is solved.

[0094] According to a second aspect of the present invention, the present invention provides a three-dimensional transient prediction device for thermal conductivity of deep formations, comprising the following modules:

[0095] The transient finite element numerical simulation module is used to obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly in the study area and set boundary conditions, and conduct transient finite element temperature numerical simulation in combination with the distribution of thermophysical parameters in the study area and the boundary conditions to obtain the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively;

[0096] The model correction amount calculation module is used to calculate the temperature field d of the underground space according to the temperature field d obs, using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm;

[0097] A model parameter updating module, used to obtain an optimal prediction model search step length satisfying the Wolfe criterion according to the prediction model correction amount Δm, and to update the prediction model parameters;

[0098] The thermal conductivity prediction module is used to carry out transient finite element temperature numerical simulation according to the updated prediction model parameters to obtain the underground space temperature field under the model parameters, calculate the data fitting difference between the measured data, and if it is less than the preset fitting difference condition, output the current prediction model parameters as the optimal thermal conductivity prediction result; otherwise, update the underground space temperature field d(m) according to the model parameters and solve the new prediction model correction value Δm.

[0099] According to a third aspect of the present invention, the present invention provides an electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the steps of the three-dimensional transient prediction method for thermal conductivity of deep formations are implemented.

[0100] The technical solution provided by the present invention has the following beneficial effects:

[0101] The present invention proposes a three-dimensional transient prediction method, device and electronic equipment for thermal conductivity of deep formations. Combining thermophysical parameters and groundwater migration laws, the Jacobian-free Krylov subspace technology based on the time chain derivative law realizes accurate and rapid prediction of thermal conductivity of deep media in the time domain. The method has strong practicality, a wide and deep prediction range, and the accuracy of the results reaches 88.97%-94.19%. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] The present invention will be further described below with reference to the accompanying drawings and embodiments, in which:

[0103] Figure 1 It is a flow chart of a three-dimensional transient prediction method of thermal conductivity of deep formations according to the present invention;

[0104] Figure 2 It is the groundwater migration rate diagram of the study area of ​​the present invention;

[0105] Figure 3 is the distribution diagram of the real thermal conductivity model of the present invention, wherein Figure 3 (a) is the thermal conductivity distribution diagram of the XZ direction cross section (y = 2.5km), Figure 3 (b) is the thermal conductivity distribution diagram of the XY direction cross section (z = 1.1 km), Figure 3 (c) is the thermal conductivity distribution diagram of the XY direction cross section (z = 3.0 km);

[0106] Figure 4 is a heat source parameter distribution diagram of the present invention;

[0107] Figure 5 It is the temperature field characteristic diagram of the underground space (y=2.5km) in the research area of ​​the present invention;

[0108] Figure 6 It is the temperature field characteristic diagram of the present invention under the prior information (y=2.5km);

[0109] Figure 7 It is the temperature field characteristic diagram of the present invention under the optimal prediction result (y=2.5km);

[0110] Figure 8 is the optimal prediction result diagram of underground thermal conductivity in the research area of ​​the present invention, where Figure 8 (a) is the thermal conductivity distribution diagram of the XZ direction cross section (y = 2.5km), Figure 8 (b) is the thermal conductivity distribution diagram of the XY direction cross section (z = 1.1 km), Figure 8 (c) is the thermal conductivity distribution diagram of the XY direction cross section (z = 3.0 km);

[0111] Fig. 9 It is a curve diagram showing the variation of the data fitting error with the number of iterations of the present invention;

[0112] Fig.10 The present invention is a structural diagram of a three-dimensional transient prediction device for thermal conductivity of deep formations. DETAILED DESCRIPTION

[0113] In order to have a clearer understanding of the technical features, purposes and effects of the present invention, specific embodiments of the present invention are now described in detail with reference to the accompanying drawings.

[0114] refer to Figure 1 , Figure 1 The present invention provides a flow chart of a method for three-dimensional transient prediction of thermal conductivity of deep formations; the method comprises the following steps:

[0115] S1: Obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly in the study area and give boundary conditions, conduct transient finite element temperature numerical simulation based on the distribution of thermophysical parameters in the study area and the boundary conditions, and obtain the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively;

[0116] In this embodiment, step S1 is specifically as follows: setting the groundwater movement rate in the x, y, and z directions of the study area (reference Figure 2 ), fluid thermal conductivity k f =0.622, the fluid migrates and immerses into the surrounding medium, causing the thermal conductivity of the medium to change, satisfying:

[0117] k=(1-φ)·k s +φ·k f

[0118] Among them, k, k s , k f are the true thermal conductivity, medium thermal conductivity and fluid thermal conductivity, respectively, with the unit of W / m·K; φ is the medium porosity, dimensionless.

[0119] A thermal conductivity anomaly is constructed in the above-mentioned uniform half-space research area, and the real thermal conductivity model distribution is as follows: Figure 3 As shown, Figure 3 (a) is the thermal conductivity distribution diagram of the XZ direction cross section (y = 2.5km), Figure 3 (b) is the thermal conductivity distribution diagram of the XY direction cross section (z = 1.1 km), Figure 3 (c) is the thermal conductivity distribution diagram of the XY direction cross section (z = 3.0 km);

[0120] refer to Figure 4 , given the heat source parameter distribution, the fluid heat capacity μ f =4.2×10 6 J / (m 3 ·K), medium heat capacity μ s =2×10 6 J / (m 3 K), and given boundary conditions:

[0121]

[0122] Among them, α=1, β=15℃ is the upper boundary condition; α=0, β=0.05W / m 2When α=0, β=0, it is the lower boundary condition; when α=0, β=0, it is the surrounding boundary condition. T is the temperature in °C, n is the direction of the boundary normal, dimensionless, and k is the thermal conductivity in W / m·K.

[0123] Under the above-mentioned thermophysical parameter distribution and boundary conditions, transient finite temperature numerical simulation was carried out to obtain the -5 The underground space temperature field after Ma is T(x,y,z)=d obs ;

[0124] As an optional implementation, transient finite element temperature numerical simulation adopts a hexahedral structure discretization method, and the temperature T, thermal conductivity k, and fluid heat capacity μ in the hexahedral unit are f 、Medium heat capacity μ s , fluid migration rate v, and heat source Q are all linearly interpolated:

[0125]

[0126]

[0127] Where i is the node number, T i , k i , μ fi , μ si 、v i , Q i They are the temperature, thermal conductivity, fluid heat capacity, medium heat capacity, fluid flow rate, and heat source of each node in the unit, respectively, in units of ℃, W / m·K, J / (m 3 ·K)、J / (m 3 ·K), m / s, W / m 3 ; N i is a shape function that satisfies:

[0128]

[0129] Among them, ξ i , η i , i is the coordinate of node i on the parent unit, dimensionless; the relationship formula between ξ, η, ζ in the parent unit and the coordinates x, y, z in the child unit is:

[0130]

[0131] Among them, x 0 ,y 0 、z 0 are the midpoints of the subunit in the x, y, and z directions respectively; a, b, and c are the side lengths of the subunit in the x, y, and z directions respectively, in m.

[0132] Given the transient heat transfer underground space temperature field differential equation is:

[0133]

[0134] Where k is the thermal conductivity, unit W / m·K, T is the temperature, unit ℃, μ f is the heat capacity of the fluid, in J / (m 3 ·K), μ s is the heat capacity of the medium, in J / (m 3 K), v is the fluid velocity, unit is m / s, Q is the heat source, unit is W / m 3 , t is time, unit is s, is the divergence, is the gradient.

[0135] The differential equation of the underground space temperature field is multiplied by the impulse function δT and integrated:

[0136]

[0137] Among them, Ω is the study area, Γ is the boundary, and δT is the impulse function.

[0138] According to the Hamiltonian operator operation rules and Ogaard formula in field theory, and substituting the boundary conditions, the first term in the above formula is rewritten as:

[0139]

[0140] Therefore, the integral equation of the temperature field of the underground space with transient heat transfer is:

[0141]

[0142] The integral equation of the underground space temperature field is discretized and linearly interpolated to obtain the unit integral of each term in the integral equation:

[0143]

[0144]

[0145]

[0146]

[0147]

[0148]

[0149] Where t is time, unit is s, δT e is the impulse function array in the unit, T eis the temperature array in the unit, K 1e ,K 2e ,K 3e ,G e are the stiffness matrices in the unit, P 1e ,P 2e is the source vector in the element, and the stiffness matrix and the coefficients in the source vector are obtained by solving:

[0150]

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157] Among them, k 1ij ,k 2ij ,k 3ij ,g ij ,p 1ij ,p 2ij The stiffness matrix K is 1e ,K 2e ,K 3e ,G e , source vector P 1e ,P 2e The elements in the , i, j, l are node labels, a, b, c are the side lengths of the subunit in the x, y, z directions, respectively, in m, ξ i , η i , i is the coordinate of node i on the parent unit.

[0158] The numbers of each node of each unit are extended to the element positions corresponding to the study area, and the linear equations of the synthesized overall stiffness matrix and source vector are obtained:

[0159]

[0160] Among them, δT is the impulse function array of the study area, T is the temperature array in the study area, K 1 ,K 2 ,K 3 , G is the stiffness matrix in the study area, P 1 ,P 2is the source vector in the study area, ignoring δT T , after simplification:

[0161]

[0162] Where K = K 1 +K 2 +K 3 , P=P 1 +P 2 .

[0163] The temperature field of the transient heat transfer underground space is differentiated at time t. After Δt, the temperature field is changed from T 1 Change to T 2 ,but In addition, the interpolated temperature within the Δt time Substituting T in the above formula, the transient temperature field linear equations are:

[0164]

[0165] By solving the above equation and realizing finite element numerical simulation, we can get the underground space temperature field T(x,y,z)=d at a given time. obs (refer to Figure 5 ), x, y, z represent the directions of x, y, and z axes respectively.

[0166] S2: According to the underground space temperature field d obs , using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton optimization algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm;

[0167] In this embodiment, step S2 is specifically as follows: setting the prior information of the research area to be a uniform half space (m ref =3.0W / m·K), and the transient finite element temperature numerical simulation was performed by substituting the distribution of thermophysical parameters and boundary conditions to obtain the underground space temperature field in the study area under the prior information (reference Figure 6 );

[0168] The objective function of Tikhonov regularization for transient temperature forward modeling is:

[0169] Φ=||D[d(m)-d obs ]||2 +λ||W(mm apr )|| 2

[0170] Where D is the data weighting matrix; d obs is the observed data, i.e., the underground space temperature field T(x, y, z); m is the current prediction model parameter vector; d(m) is the simulated data under the current prediction model parameters; W is the smooth constraint matrix; λ is the damping factor; m apr Prior information for the prediction model;

[0171] Based on the linear equations of the overall stiffness matrix and source vector in the finite element, the thermal conductivity k in the linear equations is derived as follows:

[0172]

[0173]

[0174]

[0175] Among them, K, G are the stiffness matrices in the study area, P 1 ,P 2 is the source vector in the study area, Δt is the time interval in seconds, k is the thermal conductivity array in W / m·K, T i , T i+1 are the temperature arrays at time i and i+1 respectively. At the initial time (t=0), the above formula is rewritten as:

[0176]

[0177] Among them, T 0 is the temperature array at the initial moment. To ensure that the temperature value and thermal conductivity are non-negative, the data parameters and model parameters are set as the logarithm of temperature and the logarithm of thermal conductivity during the prediction process: d = ln T, m = ln k; according to the chain rule, the elements in the Jacobian matrix are:

[0178]

[0179] Where i = 1, 2, ..., N, N is the number of data parameters, j = 1, 2, ..., M, M is the number of model parameters, d i is the data parameter of the ith node in the study area, m j represents the model parameters of the jth node in the study area, T i represents the temperature of the ith node in the study area, k j represents the thermal conductivity of the jth node in the study area, then the Jacobian matrix is ​​written as:

[0180] J=RGR'

[0181] Among them, R and R' are diagonal matrices, R'=diag(k 1 ,k 2 ,…,k M ), G is a matrix composed of elements G(i,j);

[0182] For any vector consisting of dimensionless elements v = (v 1 ,v 2 ,…,v M ) T , and the conversion formula based on the finite element stiffness matrix is:

[0183]

[0184] Among them, K(v) is the stiffness matrix when the thermal conductivity array takes v; and when t = 0, there is

[0185]

[0186] Among them, e 1 ,e 2 ,...,e M are all unit vectors, so the product of the Jacobian matrix J and any vector v can be found.

[0187] Similarly, the transpose of the Jacobian matrix is ​​calculated and the arbitrary vector y=(y 1 ,y 2 ,…,y N ) T The product of

[0188] According to the Gauss-Newton algorithm, ignoring the second-order information term, the following linear equations are obtained to solve the model correction value Δm:

[0189] HΔm=-g

[0190] Where H is the Hessian matrix, H = (J T D T DJ+λW T W), g is the gradient vector, g = [J T D T D(dd obs )+λW T W(mm apr )].

[0191] S3: according to the prediction model correction amount Δm, obtaining the optimal prediction model search step length that satisfies the Wolfe criterion, and updating the prediction model parameters;

[0192] In this embodiment, step S3 is specifically as follows: the Wolfe criterion needs to satisfy the following inequality:

[0193]

[0194] Where α is the search step length, dimensionless, and c is a dimensionless number. In the process of obtaining the search step length, let α be 1 as the initial value. If the inequality is not satisfied, let

[0195] α=ρα

[0196] Among them, ρ is a decimal greater than 0 and less than 1, and the cycle is repeated until the Wolfe criterion is met; therefore, the update formula of the model parameters is:

[0197] m k+1 =m k +αΔm k

[0198] Among them, m k 、m k+1 are the model parameters at the kth and k+1th iterations, Δm k is the model correction amount obtained in the kth iteration.

[0199] S4: According to the updated prediction model parameters, a transient finite element temperature numerical simulation is carried out to obtain the underground space temperature field under the model parameters, and the data fitting difference between the measured data and the actual data is calculated. If it is less than the preset fitting difference condition, the current prediction model parameters are the optimal thermal conductivity prediction results, and the temperature field distribution under the current model parameters is calculated as the temperature field under the optimal prediction result (reference Figure 7 ); otherwise, update the underground space temperature field d(m) according to the model parameters and solve the new prediction model correction Δm.

[0200] In this embodiment, step S4 is specifically as follows: the calculation formula of the data fitting error is as follows:

[0201]

[0202] RMS is the root mean square, which is used to measure the deviation between the measured value and the predicted value. The closer it is to 0, the smaller the deviation between the measured value and the predicted value, indicating that the prediction result is more reliable; d(m) i d obs i are the i-th predicted value and measured value respectively, n is the number of observed data, dimensionless.

[0203] The preset fitting error is 0.07. If the data fitting error is less than the preset fitting error, the iteration is stopped and the current prediction model parameters are output as the optimal prediction result (refer to Figure 8 ); otherwise, update the underground space temperature field d(m) according to the model parameters and solve the new prediction model correction Δm. Fig. 9 , Fig. 9 The curve diagram of the data fitting error of the present invention with the number of prediction iterations; in this embodiment, the RMS value between the measured value and the predicted value is calculated during each prediction process. The experimental results show that the RMS between the measured value and the predicted value converges from 1.1930 to 0.0693, and the fitting degree is good.

[0204] refer to Figure 3 , Figure 8 ,in Figure 3 is the thermal conductivity distribution diagram of the real model of the present invention, Figure 8 It is the thermal conductivity distribution diagram of the optimal prediction model of the present invention; In this embodiment, the abnormal volume overlap ratio (AVOR) is defined to evaluate the accuracy of the prediction result:

[0205]

[0206] Among them, V o is the volume of the intersection between the predicted anomaly and the real anomaly, m 3 ; V p is the volume of the predicted abnormal body, m 3 ; V t To predict abnormal bodies, m 3 The experimental results show that the AVOR of the upper thin layer anomaly reaches 88.97%, the AVOR of the lower left anomaly reaches 94.19%, and the AVOR of the lower right anomaly reaches 93.68%. This shows that the optimal prediction model has a high degree of fit with the true model.

[0207] refer to Fig.10 , Fig.10 The structure diagram of a three-dimensional transient prediction device for thermal conductivity of deep formations according to the present invention is shown in FIG. The device includes the following modules:

[0208] Transient finite element numerical simulation module 1 is used to obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly body in the study area and set boundary conditions, and carry out transient finite element temperature numerical simulation in combination with the distribution of thermophysical parameters in the study area and the boundary conditions to obtain the underground space temperature field T(x, y, z) = d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively;

[0209] Model correction amount calculation module 2, used to calculate the temperature field d of the underground space according to the temperature field d obs , using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm;

[0210] A model parameter updating module 3 is used to obtain the optimal prediction model search step length that meets the Wolfe criterion according to the prediction model correction amount Δm, and to update the prediction model parameters;

[0211] The thermal conductivity prediction module 4 is used to carry out transient finite element temperature numerical simulation according to the updated prediction model parameters to obtain the underground space temperature field under the model parameters, calculate the data fitting difference between the measured data, and if it is less than the preset fitting difference condition, output the current prediction model parameters as the optimal thermal conductivity prediction result; otherwise, update the underground space temperature field d(m) according to the model parameters and solve the new prediction model correction value Δm.

[0212] In addition, as a preferred embodiment, the present invention provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, the steps of the three-dimensional transient prediction method for thermal conductivity of deep formations are implemented, for example: obtaining the groundwater migration law and its response mechanism to thermal conductivity in a uniform half-space study area, constructing a thermal conductivity anomaly body in the study area and setting boundary conditions, and conducting transient finite element temperature numerical simulation in combination with the distribution of thermophysical property parameters in the study area and the boundary conditions to obtain the underground space temperature field T(x, y, z)=d at a given time. obs , x, y, z represent the x, y, z axis directions respectively; according to the underground space temperature field d obs, the objective function of transient temperature forward modeling is constructed by Tikhonov regularization, the Gauss-Newton algorithm is used to obtain the linear equations composed of the Hessian matrix and the gradient vector, and the Jacobian-free Krylov subspace technology based on the time chain derivation rule is used to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and the prediction model correction Δm is obtained by solving; according to the prediction model correction Δm, the optimal prediction model search step that meets the Wolfe criterion is obtained, and the prediction model parameters are updated; according to the updated prediction model parameters, the transient finite element temperature numerical simulation is carried out to obtain the underground space temperature field under the model parameters, and the data fitting difference between the calculated and measured data is calculated. If it is less than the preset fitting difference condition, the current prediction model parameters are output as the optimal thermal conductivity prediction result; otherwise, the underground space temperature field d(m) is updated according to the model parameters and the new prediction model correction Δm is solved.

[0213] It should be noted that, in this article, the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article or system including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or system. In the absence of further restrictions, an element defined by the sentence "comprises a ..." does not exclude the existence of other identical elements in the process, method, article or system including the element.

[0214] The serial numbers of the embodiments of the present invention are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments. In a unit claim that lists several means, several of these means may be embodied by the same hardware item. The use of the words first, second, and third, etc. does not indicate any order and these words may be interpreted as identifiers.

[0215] The above are only preferred embodiments of the present invention, and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification and drawings, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A three-dimensional transient prediction method for thermal conductivity of deep formations, characterized in that: The following steps are involved: Obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly in the study area and set boundary conditions, conduct transient finite element temperature numerical simulation based on the distribution of thermophysical parameters in the study area and the boundary conditions, and obtain the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively; According to the underground space temperature field d obs , using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm; According to the prediction model correction amount Δm, an optimal prediction model search step length satisfying the Wolfe criterion is obtained, and prediction model parameters are updated; According to the updated prediction model parameters, transient finite element temperature numerical simulation is carried out to obtain the underground space temperature field under the model parameters, and the data fitting difference between the measured data is calculated. If it is less than the preset fitting difference condition, the current prediction model parameters are output as the optimal thermal conductivity prediction result; otherwise, the underground space temperature field d(m) is updated according to the model parameters and the new prediction model correction Δm is solved.

2. The three-dimensional transient prediction method for thermal conductivity of deep formations according to claim 1, characterized in that: The groundwater migration law and its response mechanism to thermal conductivity include: Assume that the groundwater movement rates in the x, y, and z directions of the study area are v x 、v y 、v z , the fluid migrates and immerses into the surrounding medium, causing the thermal conductivity of the medium to change, satisfying: k=(1-φ)·k s +φ·k f Among them, k, k s , k f are the true thermal conductivity, medium thermal conductivity and fluid thermal conductivity, respectively, with the unit of W / m·K; φ is the medium porosity, dimensionless.

3. The three-dimensional transient prediction method for thermal conductivity of deep formations according to claim 1, characterized in that: The boundary conditions include: Among them, α and β are constants, α=1, β=T up When is the upper boundary condition; α=0, β=q down When α=0, β=0, it is the lower boundary condition; when α=0, β=0, it is the surrounding boundary condition; T is the temperature, unit is ℃, T up is the surface temperature, unit is °C, n is the direction of the boundary normal, dimensionless, k is the thermal conductivity, unit is W / m·K, q down is the heat flow value, unit is W / m 2 .

4. The three-dimensional transient prediction method for thermal conductivity of deep formations according to claim 1, characterized in that: The transient finite element temperature numerical simulation adopts a hexahedral structure discretization method, and the temperature T, thermal conductivity k, and fluid heat capacity μ in the hexahedral unit are f 、Medium heat capacity μ s , fluid migration rate v, and heat source Q are all linearly interpolated: Where i is the node number, T i , k i , μ fi , μ si 、v i , Q i They are the temperature, thermal conductivity, fluid heat capacity, medium heat capacity, fluid flow rate, and heat source of each node in the unit, respectively, in units of ℃, W / m·K, J / (m 3 ·K)、J / (m 3 ·K), m / s, W / m 3 ; N i is a shape function that satisfies: Among them, ξ i , η i , i is the coordinate of node i on the parent unit, dimensionless; the relationship formula between ξ, η, ζ in the parent unit and the coordinates x, y, z in the child unit is: Among them, x0, y0, z0 are the midpoints of the subunit in the x, y, z directions respectively; a, b, c are the side lengths of the subunit in the x, y, z directions respectively, in m.

5. The three-dimensional transient prediction method of thermal conductivity of deep formations according to claim 1, characterized in that: The step of carrying out transient finite element temperature numerical simulation in combination with the distribution of thermophysical parameters in the study area and the boundary conditions to obtain the temperature field of the underground space at a given time includes: Given the transient heat transfer underground space temperature field differential equation is: Where k is the thermal conductivity, unit W / m·K, T is the temperature, unit ℃, μ f is the heat capacity of the fluid, in J / (m 3 ·K), μ s is the heat capacity of the medium, in J / (m 3 K), v is the fluid velocity, unit is m / s, Q is the heat source, unit is W / m 3 , t is time, unit is s, is the divergence, is the gradient; The differential equation of the underground space temperature field is multiplied by the impulse function δT and integrated: Among them, Ω is the study area, Γ is the boundary, and δT is the impulse function; According to the Hamiltonian operator operation rules and Ogaard formula in field theory, and substituting the boundary conditions, the first term in the above formula is rewritten as: The integral equation of the temperature field of the underground space under transient heat transfer is obtained as follows: The integral equation of the temperature field of the underground space of transient heat transfer is discretized and linearly interpolated to obtain the unit integral of each term in the integral equation: Where t is time, unit is s, δT e is the impulse function array in the unit, T e is the temperature array in the unit, K 1e ,K 2e ,K 3e ,G e is the stiffness matrix within the element, P 1e ,P 2e is the source vector in the element, and the stiffness matrix and the coefficients in the source vector are obtained by solving: Among them, k 1ij ,k 2ij ,k 3ij ,g ij ,p 1ij ,p 2ij The stiffness matrix K is 1e ,K 2e ,K 3e ,G e and the source vector P 1e ,P 2e The elements in the , i, j, l are node labels, a, b, c are the side lengths of the subunit in the x, y, z directions, respectively, in m, ξ i , η i , i is the coordinate of node i on the parent unit; The numbers of each node of each unit are extended to the element positions corresponding to the study area, and the linear equations of the synthesized overall stiffness matrix and source vector are obtained: Among them, δT is the impulse function array of the study area, T is the temperature array in the study area, K1, K2, K3, G are the stiffness matrices in the study area, P1, P2 are the source vectors in the study area, and δT is omitted. T , after simplification: Among them, K = K1 + K2 + K3, P = P1 + P2; The temperature field of the transient heat transfer underground space is differentiated at time t. After Δt, the temperature field changes from T1 to T2. Then In addition, the interpolated temperature within the Δt time Substituting T in the above formula, the transient temperature field linear equations are: By solving the above equation and realizing finite element numerical simulation, we can get the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively.

6. The three-dimensional transient prediction method of thermal conductivity of deep formations according to claim 5, characterized in that: The objective function of the Tikhonov regularization to construct the transient temperature forward modeling is: Φ=||D[d(m)-d obs ]|| 2 +λ||W(m-m apr )|| 2 Where D is the data weighting matrix; d obs is the observed data, i.e., the underground space temperature field T(x, y, z); m is the current prediction model parameter vector; d(m) is the simulated data under the current prediction model parameters; W is the smooth constraint matrix; λ is the damping factor; m apr Prior information for the prediction model; Based on the linear equations of the overall stiffness matrix and source vector in the finite element, the thermal conductivity k in the linear equations is derived as follows: Among them, K, G are the stiffness matrices in the study area, P1, P2 are the source vectors in the study area, Δt is the time interval in seconds, k is the thermal conductivity array in W / m·K, T i 、T i+1 are the temperature arrays at time i and i+1 respectively; at the initial time (t=0), the above formula is rewritten as: Among them, T0 is the temperature array at the initial moment; in the prediction process, the data parameters and model parameters are set as the logarithm of temperature and the logarithm of thermal conductivity: d = lnT, m = lnk; according to the chain rule, the elements in the Jacobian matrix are: Where i = 1, 2, ..., N, N is the number of data parameters, j = 1, 2, ..., M, M is the number of model parameters, d i is the data parameter of the ith node in the study area, m j represents the model parameters of the jth node in the study area, T i represents the temperature of the ith node in the study area, k j represents the thermal conductivity of the jth node in the study area, then the Jacobian matrix is ​​written as: J=RGR' Among them, R and R' are diagonal matrices, R'=diag(k1,k2,…,k M ), G is a matrix composed of elements G(i,j); For any vector v consisting of dimensionless elements, M ) T , and the conversion formula based on the finite element stiffness matrix is: Among them, K(v) is the stiffness matrix when the thermal conductivity array takes v; and when t = 0, there is Among them, e1,e2,...,e M All are unit vectors, find the product of the Jacobian matrix J and any vector v; Similarly, the transpose of the Jacobian matrix and any vector y=(y1,y2,…,y N ) T The product of According to the Gauss-Newton algorithm, ignoring the second-order information term, the following linear equations are obtained to solve the model correction value Δm: HΔm=-g Where H is the Hessian matrix, H = (J T D T DJ+λW T W), g is the gradient vector, g = [J T D T D(dd obs )+λW T W(mm apr )].

7. The three-dimensional transient prediction method of thermal conductivity of deep formations according to claim 6, characterized in that: The optimal prediction model search step length satisfying the Wolfe criterion is obtained according to the prediction model correction amount Δm, which needs to satisfy the following inequality: Where α is the search step length, dimensionless, and c is a dimensionless number. In the process of obtaining the search step length, let α be 1 as the initial value. If the inequality is not satisfied, let α=ρα Among them, ρ is a decimal greater than 0 and less than 1, and the cycle is repeated until the Wolfe criterion is met; the update formula of the model parameters is: m k+1 =m k +αΔm k Among them, m k 、m k+1 are the model parameters at the kth and k+1th iterations, Δm k is the model correction amount obtained in the kth iteration.

8. The three-dimensional transient prediction method for thermal conductivity of deep formations according to claim 1, characterized in that: The calculation formula of the data fitting error is as follows: RMS is the root mean square, which is used to measure the deviation between the measured value and the predicted value; d(m) i ,d obs i are the ith predicted value and measured value respectively, n is the number of observed data, dimensionless; If the data fitting error is less than the preset fitting error, the iteration is stopped and the current prediction model parameters are output as the optimal prediction result; otherwise, the underground space temperature field d(m) is updated according to the model parameters and the new prediction model correction Δm is solved.

9. A three-dimensional transient prediction device for thermal conductivity of deep formations, characterized in that: Includes the following modules: The transient finite element numerical simulation module is used to obtain the groundwater migration law and its response mechanism to thermal conductivity in the uniform half-space study area, construct a thermal conductivity anomaly in the study area and set boundary conditions, and conduct transient finite element temperature numerical simulation in combination with the distribution of thermophysical parameters in the study area and the boundary conditions to obtain the underground space temperature field T(x,y,z)=d at a given time. obs , x, y, z represent the directions of x, y, and z axes respectively; The model correction amount calculation module is used to calculate the temperature field d of the underground space according to the temperature field d obs , using Tikhonov regularization to construct the objective function of transient temperature forward modeling, using Gauss-Newton algorithm to obtain the linear equations consisting of Hessian matrix and gradient vector, using Jacobian-free Krylov subspace technology based on time chain derivation law to solve the Jacobian matrix of the first-order derivative of the temperature numerical simulation response to the prediction model parameters in the Hessian matrix and the product of the transpose of the Jacobian matrix and an arbitrary vector, and solve to obtain the prediction model correction Δm; A model parameter updating module, used to obtain an optimal prediction model search step length satisfying the Wolfe criterion according to the prediction model correction amount Δm, and to update the prediction model parameters; The thermal conductivity prediction module is used to carry out transient finite element temperature numerical simulation according to the updated prediction model parameters to obtain the underground space temperature field under the model parameters, calculate the data fitting difference between the measured data, and if it is less than the preset fitting difference condition, output the current prediction model parameters as the optimal thermal conductivity prediction result; otherwise, update the underground space temperature field d(m) according to the model parameters and solve the new prediction model correction value Δm.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the steps of the three-dimensional transient prediction method of thermal conductivity of deep formations as described in any one of claims 1-8 are implemented.

Images