Physical auxiliary machine learning deep geothermal heat exchanger operating state prediction method

Abstract

The application relates to the technical field of deep geothermal energy development and utilization, in particular to a physical auxiliary machine learning deep geothermal heat exchanger operation state prediction method, which comprises the following steps: acquiring wellbore structure parameters, drilling exploration data and operation condition parameters; an analytical heat transfer model considering rock-soil layering and a ground temperature gradient is established according to the drilling exploration data, and an initial ground temperature field is determined; the integral path of a finite line heat source is subjected to layering treatment, equivalent thermal physical property parameters are calculated, and a layering heat response relationship is established; a heat response matrix is calculated based on the layering heat response relationship, the initial ground temperature field and a historical heat flow increment are combined to solve well wall temperature, and an analytical heat transfer model output is obtained; simulation operation data are generated by using the analytical heat transfer model, and a simulation data set is constructed; a physical auxiliary machine learning model is trained by using the simulation data set and measured operation data; the operation condition parameters to be predicted and the analytical heat transfer model output are input into the trained physical auxiliary machine learning model, and an operation state prediction result is obtained.

Description

Technical Field

[0001] This application relates to the field of deep geothermal energy development and utilization technology, and in particular to a physical-assisted machine learning method for predicting the operating status of deep geothermal heat exchangers. Background Technology

[0002] With the continuous development of geothermal energy development and utilization technology, deep geothermal heat exchangers, due to their ability to utilize the relatively high and stable underground geothermal resources for building heating, cooling, and integrated energy utilization, have been widely applied in building energy systems, district energy supply, and industrial waste heat utilization. During the operation of deep geothermal heat exchangers, continuous heat exchange occurs between the well fluid and the surrounding rock. This heat transfer process is simultaneously influenced by factors such as well structure, fluid operating conditions, soil and rock stratification characteristics, differences in formation thermal properties, and geothermal gradient distribution. For deep geothermal heat exchangers with significant well depth and heterogeneous formations, the underground heat transfer process typically exhibits obvious multi-scale, strongly coupled, and long-period characteristics. Therefore, accurate prediction of the operating status of deep geothermal heat exchangers is a crucial foundation for achieving optimized system operation, performance evaluation, and lifespan management.

[0003] Currently, the prediction of the operating status of deep geothermal heat exchangers typically employs analytical models, equivalent thermal resistance models, three-dimensional numerical simulation methods, and data-driven prediction methods. Among these, modeling deep geothermal heat exchangers involves multi-media coupled heat transfer processes involving underground rock and soil, wellbore fluids, casing, central pipe, and backfill materials. These wells are typically thousands of meters deep, and their operating timescales often span several years. While traditional high-fidelity three-dimensional numerical simulation methods can accurately describe complex underground heat transfer processes, they suffer from large computational scales, high computational costs, and long computation cycles, making them unsuitable for applications requiring extensive parameter combination analysis, scheme optimization design, and long-term operational prediction. Analytical models, while offering advantages such as fast computation speed and ease of engineering application, are mostly based on assumptions such as homogeneous strata, constant thermophysical parameters, or ideal boundary conditions. They fail to adequately consider complex heat transfer processes such as rock and soil stratification, geothermal gradient changes, and long-term heat accumulation, making it difficult to accurately reflect the actual heat transfer characteristics under complex deep geological conditions. For example, patent CN119416596A provides a method for predicting the heat extraction capacity of coaxial casing heat exchangers in deep geothermal wells; patent CN120805781A provides a reservoir-wellbore integrated prediction method for geothermal extraction; and patent CN120951705A provides a wellbore thermo-coupling simulation method, device, and equipment based on a fractional-order model.

[0004] On the other hand, in recent years, there has been a growing trend of using machine learning to build predictive models for the operating status of deep geothermal heat exchangers. These methods utilize historical operating data to train the predictive model, improving the speed of operating status prediction. However, purely data-driven models typically rely on a large amount of well-covered and evenly distributed measured operating data. Deep geothermal heat exchangers have long testing cycles and high operating costs, resulting in a limited amount of available measured operating data, which is difficult to cover different wellbore structures, formation conditions, and operating conditions. Therefore, under small sample conditions, existing data-driven models are prone to problems such as insufficient generalization ability, weak extrapolation ability, and poor adaptability to complex formation conditions. Furthermore, due to the lack of constraints on underground heat transfer mechanisms, their prediction results are difficult to maintain physical consistency, easily leading to significant prediction biases in long-term operating predictions and predictions of unknown operating conditions, making it difficult to balance prediction accuracy, computational efficiency, and engineering applicability. Summary of the Invention

[0005] (a) Technical problems to be solved In view of the above-mentioned shortcomings and deficiencies of the prior art, this application provides a physical-aided machine learning method for predicting the operating status of deep geothermal heat exchangers. It solves the technical problems of insufficient prediction accuracy of analytical heat transfer models under complex stratified geological conditions, low computational efficiency of traditional numerical simulation, high dependence of data-driven models on measured samples, and insufficient accuracy and generalization ability of operating status prediction under complex operating conditions in the prior art.

[0006] (II) Technical Solution To achieve the above objectives, the main technical solutions adopted in this application include: This application provides a method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning, comprising: S1, acquiring wellbore structural parameters, drilling exploration data, and operating condition parameters of the deep geothermal heat exchanger; establishing an analytical heat transfer model considering soil-rock stratification and geothermal gradient based on the drilling exploration data, and determining the initial geothermal field of the analytical heat transfer model based on the geothermal gradient of each stratum; S2, performing stratification processing on the finite linear heat source integral path based on the analytical heat transfer model, and calculating the corresponding equivalent thermal property parameters based on the thermal property parameters of each stratum, thereby establishing a stratified thermal response relationship that satisfies the interlayer temperature continuity condition and the heat flow continuity condition; S3, calculating the thermal response matrix based on the stratified thermal response relationship, and solving for the wellbore temperature by combining the initial geothermal field and historical heat flow increments; establishing the coupled heat transfer relationship between the wellbore fluid and the surrounding rock based on the wellbore temperature, and obtaining the analytical heat transfer... Model output; S4. Using the analytical heat transfer model, perform parametric calculations under different wellbore structural parameters, operating condition parameters, and drilling exploration data conditions to generate simulated operating data. Based on the simulated operating data, perform adaptive enhancement to construct a simulated dataset; S5. Use the simulated dataset and measured operating data to train a physical-aided machine learning model and obtain the trained physical-aided machine learning model; S6. Obtain the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted. Based on the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted, obtain the output of the analytical heat transfer model corresponding to the operating condition to be predicted. Input the wellbore structural parameters, drilling exploration data, operating condition parameters, and the output of the analytical heat transfer model into the trained physical-aided machine learning model to obtain the predicted operating status of the deep geothermal heat exchanger.

[0007] (III) Beneficial Effects The physical-aided machine learning method for predicting the operating status of deep geothermal heat exchangers provided in this application establishes an analytical heat transfer model that considers soil and rock stratification and geothermal gradients. Combined with the calculation of the stratified thermal response relationship and thermal response matrix of a finite linear heat source, it achieves rapid analytical calculation of the coupled heat transfer process between wellbore fluid and surrounding rock under complex stratified geological conditions, thus improving the computational accuracy of the analytical heat transfer model. Furthermore, by generating simulated operating data using the analytical heat transfer model, adaptively enhancing the simulated data, and training the physical-aided machine learning model with measured operating data, it achieves a synergistic integration of physical mechanisms and data-driven approaches. This reduces the dependence on a large number of measured samples while improving the model's prediction accuracy, generalization ability, and computational efficiency under complex geological conditions and different operating conditions, thereby enabling rapid, accurate, and reliable prediction of the operating status of deep geothermal heat exchangers. Attached Figure Description

[0008] Figure 1This is a flowchart illustrating a method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning according to an embodiment of this application. Figure 2 This is a schematic diagram of the wellbore structure and layered formations suitable for a deep coaxial casing geothermal heat exchanger according to an embodiment of this application. Detailed Implementation

[0009] To better explain and facilitate understanding of this application, the following detailed description of the application is provided in conjunction with the accompanying drawings and specific embodiments.

[0010] Example 1; Figure 1 This is a flowchart illustrating a method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning, according to an embodiment of this application. Figure 1 As shown, the physical-aided machine learning method for predicting the operating status of deep geothermal heat exchangers includes: S1. Obtain the wellbore structural parameters, drilling exploration data, and operating condition parameters of the deep geothermal heat exchanger; establish an analytical heat transfer model considering soil and rock stratification and geothermal gradient based on the drilling exploration data, and determine the initial geothermal field of the analytical heat transfer model based on the geothermal gradient of each stratum. The wellbore structural parameters include at least the well depth, well diameter, well wall radius, wellbore heat exchange section length, central tube inner diameter, central tube outer diameter, annulus inner diameter, annulus outer diameter, casing inner diameter, casing outer diameter, central tube wall thickness, and casing wall thickness; the drilling exploration data includes at least the top and bottom depths of each formation, layer thickness, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradient; the operating condition parameters include at least the inlet temperature, circulation flow rate, heat load, and operating time. Wherein, the initial geothermal field satisfies: ; in, For depth The initial ground temperature at that location, This is the surface reference temperature. For the first Geothermal gradient of the strata For the first Geothermal gradient of the strata For the first The thickness of the strata, and The first The depth of the top and bottom interfaces of the strata; The initial geothermal field satisfies the temperature continuity condition at each stratum interface: ; Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time; Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time; where the superscript "-" indicates approaching the layer from the top side of the interface, and the superscript "+" indicates approaching the interface from the bottom side of the interface; And satisfy the heat flow continuity condition: ; in, Indicates the first The thermal conductivity of the strata Indicates the first The thermal conductivity of the strata, Indicates from the first One side of the layer approaches the interface Temperature gradient over time; Indicates from the first One side of the layer approaches the interface The temperature gradient over time; n is the normal direction of the interlayer interface; Indicates from the first One side of the layer approaches the interlayer interface. Indicates from the first One side of the layer approaches the interlayer interface.

[0011] For example, in this embodiment of the application, the wellbore structural parameters, drilling exploration data, and operating condition parameters of the deep geothermal heat exchanger are first obtained. The wellbore structural parameters are used to determine the wellbore geometry and fluid flow channels; the drilling exploration data is used to characterize the underground strata distribution and thermal properties; and the operating condition parameters are used to describe the actual operating state of the heat exchanger. For instance, for a deep geothermal heat exchange well with a depth of approximately 3000m, wellbore structural parameters such as well depth, well diameter, well wall radius, wellbore heat exchange section length, central tube inner diameter, central tube outer diameter, annulus inner diameter, annulus outer diameter, casing inner diameter, casing outer diameter, central tube wall thickness, and casing wall thickness can be obtained. Simultaneously, formation data formed during drilling is acquired, including drilling exploration data such as the top and bottom depths of each stratum, layer thickness, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradient. Operating condition parameters may include parameters such as inlet temperature, circulation flow rate, heat load, and continuous operating time. Furthermore, an analytical heat transfer model considering the effects of soil and rock stratification and geothermal gradients is established based on the obtained drilling exploration data. Since deep geothermal wells typically traverse multiple soil and rock layers with different properties, and the thermal conductivity and geothermal gradients of each stratum vary significantly, this application no longer adopts the homogeneous stratification assumption. Instead, based on the geological information obtained from drilling exploration, the underground surrounding rock is divided into multiple strata, and each is assigned corresponding thermal properties such as thermal conductivity, density, specific heat capacity, and thermal diffusivity, thereby constructing an analytical heat transfer model that reflects the actual underground stratification characteristics. After establishing the analytical heat transfer model, the initial underground geothermal field is determined based on the geothermal gradients corresponding to each stratum. Specifically, using the surface reference temperature as the initial boundary condition, when the depth is within the i-th stratum, the initial geothermal temperature at depth z is determined by the cumulative geothermal gradients of the overlying strata and the geothermal gradient of the current stratum, i.e.: The first term represents the surface reference temperature, the second term represents the cumulative temperature rise of the overlying strata, and the third term represents the temperature rise within the current stratum with increasing depth. Through the above calculations, the initial geothermal distribution, which varies continuously along the entire well depth, can be obtained.

[0012] For example, when drilling through three strata, with the first layer being 500m thick and having a geothermal gradient of 2.8℃ / 100m, the second layer being 1000m thick and having a geothermal gradient of 3.4℃ / 100m, and the third layer being 1500m thick and having a geothermal gradient of 4.0℃ / 100m, and the surface reference temperature being 18℃, the initial geothermal temperature at a depth of 2000m in the third layer can be calculated from the cumulative temperature rise of the first two layers and the local temperature rise of the third layer, without the need for approximate calculation using a uniform geothermal gradient. This allows for a more realistic reflection of the underground temperature field distribution. Furthermore, to ensure that the analytical heat transfer model satisfies actual heat transfer laws, the initial geothermal field established in this application also satisfies the temperature continuity condition and heat flow continuity condition at the interlayer interfaces. Specifically, at any stratum interface, the initial geothermal temperatures on both the upper and lower sides remain continuous, i.e., satisfying: This indicates that there is no abrupt change in underground temperature at the interlayer interface; at the same time, the heat flow continuity condition is satisfied on both sides of the interface. This means that the heat flux density transferred on both sides of the interface remains consistent. Different strata with varying thermal conductivity correspond to different temperature gradients, thus ensuring that interlayer heat conduction satisfies the law of energy conservation. The initial geothermal field established through these constraints not only accurately reflects the stratified geological conditions underground but also provides accurate initial boundary conditions for subsequent calculations of finite linear heat sources, construction of thermal response matrices, and wellbore fluid coupling heat transfer, thereby improving the computational accuracy of the entire analytical heat transfer model.

[0013] S2. Based on the analytical heat transfer model, the finite linear heat source integral path is processed into layers, and the corresponding equivalent thermal properties are calculated according to the thermal property parameters of each layer, thereby establishing a layered thermal response relationship that satisfies the interlayer temperature continuity condition and the heat flow continuity condition; S2 specifically includes: S21. Determine the finite linear heat source integration path based on the positional relationship between the wellbore heat exchange section and the layered formations; S22. When the finite linear heat source integration path crosses a formation interface, the finite linear heat source integration path is segmented according to the formation interface to form multiple sub-paths located in different formations. S23. Calculate the equivalent thermal properties of each sub-path based on the thermal property parameters of the formation corresponding to each sub-path; wherein, the thermal property parameters include thermal conductivity, density, and specific heat capacity; the equivalent thermal property parameters include: equivalent thermal conductivity. equivalent density Equivalent specific heat capacity and equivalent thermal diffusivity ;in, ; ; ; ; To indicate the first The proportion of each sub-path length to the integral path length of the finite line heat source. , and They represent the first The thermal conductivity, density, and specific heat capacity of the strata corresponding to each sub-path; Indicates the subpath number involved in the summation; Indicates the first The proportion of the length of each sub-path to the length of the integral path of the finite linear heat source; Indicates the first The density of the stratum where each sub-path is located; S24. Calculate the temperature response function of the finite linear heat source in the corresponding sub-path based on the equivalent thermal property parameters of each sub-path. Combine the temperature response functions according to the position order of each sub-path in the integral path of the finite linear heat source, and satisfy the temperature continuity condition and heat flow continuity condition of the adjacent strata interface to establish the layered thermal response relationship of the finite linear heat source in the layered strata.

[0014] For example, in this embodiment, after establishing an analytical heat transfer model considering soil stratification and geothermal gradient, it is necessary to further establish the heat transfer response relationship of the finite linear heat source in the stratified strata. Since deep geothermal heat exchangers are typically buried thousands of meters underground, the wellbore heat exchange section often traverses multiple different rock strata simultaneously, and the thermal conductivity, heat storage capacity, and heat diffusion capacity of different rock strata vary significantly. If the homogeneous strata assumption in the traditional finite linear heat source model is still adopted, the thermophysical parameters along the entire integration path will be considered consistent, failing to truly reflect the heat transfer differences generated when heat propagates across different rock strata, thereby reducing the accuracy of subsequent thermal response matrix calculations. Therefore, this application first performs stratification processing on the finite linear heat source integration path, and then establishes a stratified thermal response relationship that satisfies the inter-layer continuity condition. Specifically, the finite linear heat source integration path is first determined based on the positional relationship between the wellbore heat exchange section and the underground strata. The finite linear heat source integration path is the length of the heat source corresponding to the integration calculation of the finite linear heat source model, which is usually consistent with the length of the wellbore heat exchange section. For example, a deep geothermal well has an effective heat exchange section length of 2500m. Drilling exploration results indicate that this heat exchange section successively traverses sandstone, limestone, and granite layers. Therefore, the finite linear heat source integration path corresponds to a depth range of 0–2500m. It is then determined whether the finite linear heat source integration path crosses a formation interface. When the finite linear heat source integration path passes through multiple formations, each formation interface is used as a dividing point, and the original integration path is divided into multiple sub-paths located in different formations. For example, in the aforementioned 2500m heat exchange section, sandstone layers are distributed from 0–800m, limestone layers from 800–1700m, and granite layers from 1700–2500m. Therefore, the finite linear heat source integration path is divided into three sub-paths: the first sub-path is 800m long, the second sub-path is 900m long, and the third sub-path is 800m long, each sub-path corresponding to its respective formation.

[0015] For each sub-path, the equivalent thermal properties of that sub-path are calculated based on the thermal properties of the rock stratum in which it is located. These thermal properties include at least thermal conductivity, density, and specific heat capacity. For example, the parameters corresponding to the three strata mentioned above could be: First layer sandstone: thermal conductivity 2.5 W / (m·K), density 2200 kg / m³, specific heat capacity 900 J / (kg·K); Second layer limestone: thermal conductivity 3.2 W / (m·K), density 2500 kg / m³, specific heat capacity 850 J / (kg·K); Third layer granite: thermal conductivity 3.8 W / (m·K), density 2650 kg / m³, specific heat capacity 790 J / (kg·K). The length weight is calculated based on the proportion of each sub-path length to the total integral path length, i.e.: =800 / 2500=0.32; =900 / 2500=0.36; =800 / 2500=0.32; then calculate the equivalent thermal properties parameters corresponding to the entire finite linear heat source integration path. Through calculation, the equivalent thermal conductivity and thermal storage performance of the finite linear heat source integration path after comprehensively considering multiple different strata can be obtained. Further, based on the equivalent thermal properties parameters corresponding to each sub-path, the temperature response function of the finite linear heat source in the corresponding sub-path is calculated separately. Since each sub-path corresponds to a different stratum, the transient heat transfer response of each sub-path is calculated using its own corresponding equivalent thermal conductivity and equivalent thermal diffusivity, instead of using uniform stratum parameters. For example, in the sandstone layer, the equivalent thermal conductivity and equivalent thermal diffusivity corresponding to the first sub-path are used to calculate the temperature response; when the heat source propagates into the limestone layer, the calculation automatically switches to the equivalent thermal properties parameters corresponding to the second sub-path; after entering the granite layer, the parameters corresponding to the third sub-path are used for calculation. Then, according to the spatial order of each sub-path in the finite linear heat source integration path, the temperature response functions corresponding to each sub-path are sequentially connected and combined. Simultaneously, temperature continuity and heat flow continuity conditions are applied at the boundaries of various strata to ensure that the combined overall temperature response function remains physically continuous at the interlayer interfaces. Specifically, the temperature continuity condition guarantees that there are no abrupt temperature changes on both sides of the interface, while the heat flow continuity condition ensures that the heat flow density remains consistent on both sides of the interface. These two continuity conditions allow for a natural transition in temperature propagation from a finite linear heat source across different rock strata, avoiding artificial jumps in the temperature field caused by directly switching thermophysical parameters.

[0016] Finally, by sequentially combining the temperature response functions of the sub-paths and satisfying the temperature continuity and heat flow continuity conditions at each interface, the layered thermal response relationship of the finite linear heat source in the layered strata is established. Compared with the traditional homogeneous finite linear heat source model, the layered thermal response relationship established in this application can fully consider the differences in thermal properties of different underground rock strata and the interlayer heat transfer characteristics, and more realistically reflect the heat propagation law in complex layered strata during the long-term operation of deep geothermal heat exchangers. This provides a more accurate physical basis for subsequent thermal response matrix calculation, wellbore temperature solution, and wellbore fluid coupling heat transfer analysis, thereby improving the prediction accuracy of the entire operating state prediction model.

[0017] S3. Calculate the thermal response matrix based on the layered thermal response relationship, and solve the wellbore temperature by combining the initial geothermal field and historical heat flow increment; establish the coupled heat transfer relationship between the wellbore fluid and the surrounding rock based on the wellbore temperature, and obtain the output of the analytical heat transfer model; In this embodiment, after obtaining the layered thermal response relationship established in step S2, the wellbore temperature is further calculated, and the coupled heat transfer relationship between the wellbore fluid and the surrounding rock is established to obtain the analytical heat transfer model output. The thermal response relationship described in this application represents the temperature response function of different locations on the wellbore under unit heat flux from a finite linear heat source. It comprehensively considers different formation thermal conductivity, formation thickness, geothermal gradient, and interlayer temperature and heat flux continuity conditions, thus reflecting the true heat transfer law in deep heterogeneous formations.

[0018] Since the wellbore undergoes a continuous heat exchange process during actual operation, the heat injected underground at different times will continuously affect the wellbore temperature at subsequent times. Therefore, this application further utilizes the thermal response relationship to establish a thermal response matrix, and combines historical heat flow information to calculate the wellbore temperature, and then further solves for the wellbore fluid temperature distribution. Step S3 specifically includes the following steps S31 to S34: S31. Based on the hierarchical thermal response relationship, the transient heat transfer response of the finite linear heat source under different time lags is pre-calculated. The thermal response functions between each calculation node in the wellbore and each heat source segment are calculated, and a thermal response matrix is ​​constructed according to the thermal response functions corresponding to each time lag. For example, in step S31, the hierarchical thermal response relationship obtained in step S2 is first used to pre-calculate the finite linear heat source under multiple time lag conditions. Here, the time lag in this application represents the time interval experienced by a unit heat flux from the application time to the current calculation time. For example, if the time step is... Taking 1 hour as an example, the time lag can correspond to multiple historical moments such as 1 hour, 2 hours, 3 hours, ... 100 hours. Since underground heat conduction is a typical diffusion process, the temperature influence range corresponding to the same unit heat flow is different at different time lags. Therefore, it is necessary to calculate the temperature response for each time lag separately.

[0019] Specifically, the entire wellbore is discretized into multiple computational nodes along the depth direction, and the finite linear heat source integration path is discretized into multiple heat source segments. For example, a 3000m deep geothermal well can be discretized into computational nodes of 20m each, resulting in a total of 150 wellbore computational nodes; at the same time, the entire heat exchange section is also divided into 150 finite linear heat source segments.

[0020] Subsequently, for each time lag Calculate separately: the first The wellbore computing node is affected by the first The temperature change resulting from a unit heat flux in a heat source segment is called the thermal response function. The thermal response function represents the temperature increase at the corresponding computational node per unit heat flux input. Due to the existence of surface boundaries underground, this application considers both real heat sources and mirrored surface heat sources. Real heat sources refer to the heat sources formed by the actual heat released from the wellbore. Mirrored heat sources are established using the mirror method, and their purpose is to satisfy surface boundary conditions, making the underground temperature field more consistent with actual engineering conditions.

[0021] For example: well wall radius =0.09m; Calculate node depth =1800m; the integral variable h of the heat source segment is 1700m. Therefore: the actual distance to the heat source is: Distance from mirrored heat source: ;because: Significantly greater than Therefore, the influence of the mirrored heat source is relatively weak, but it can still improve the accuracy of the model boundary. Subsequently, the equivalent thermal conductivity and equivalent thermal diffusivity obtained in step S2 are used to calculate the transient temperature response corresponding to the real heat source and the mirrored heat source, respectively. Here, the equivalent thermal conductivity determines the heat transfer capability, and the equivalent thermal diffusivity determines the heat diffusion rate.

[0022] Since the heat source integration path may cross multiple strata, different integration locations correspond to different equivalent thermophysical parameters. Therefore, integration along the entire integration path is required. After integration, we obtain: the... The calculation point and the first There is a time lag between the heat source segments The corresponding thermal response matrix elements are as follows: Its physical meaning can be understood as: the first After the heat flow in the section releases 1 W / m of heat, The moment makes the first The calculation determines the temperature rise in degrees Celsius at each calculation point. This process is then repeated for all calculation points, all heat source segments, and all time lags. For example, with 150 calculation points, 150 heat source segments, and 300 time steps, 300 sets of 150×150 thermal response matrices are generated. All these matrices are pre-calculated and stored in memory for quick subsequent retrieval. Therefore, repeated integration is unnecessary during the prediction phase, significantly reducing computation time.

[0023] S32. Based on the initial geothermal field, the thermal response matrix, and the historical heat flux increment, calculate the wellbore temperature corresponding to each time step by recursive calculation through historical heat flux convolution. In step S32, the wellbore temperature is calculated using the thermal response matrix pre-established in step S31. Because the surrounding underground rock has a large heat capacity and a low heat diffusion rate, its heat transfer process exhibits significant thermal inertia. That is, the current wellbore temperature is affected not only by the heat flow input at the current time step but also by the combined influence of the continuous diffusion of heat flow from previous time steps. Therefore, this application employs a historical heat flow convolution recursive method to calculate the wellbore temperature, fully considering the time accumulation effect in the underground heat transfer process, thereby more realistically reflecting the wellbore temperature variation pattern during the long-term operation of the deep geothermal heat exchanger.

[0024] Specifically, the heat flow per unit well depth along the depth direction at the nth time step is represented as a well depth heat flow vector. (That is, the heat transfer per unit well depth corresponding to each calculation node at the nth time step is arranged according to the well depth direction to form a well depth heat flow vector) The well depth heat flow vector consists of the heat transfer per unit well depth corresponding to each depth calculation node in the wellbore. Assuming the wellbore is discretized into 150 calculation points along the depth direction, the well depth heat flow vector contains 150 components, of which... This represents the heat transfer per unit well depth corresponding to the 20th calculation node. Because the heat transfer intensity between the well fluid and the surrounding rock changes continuously with variations in inlet temperature, circulation flow rate, and underground temperature field during the operation of a deep geothermal heat exchanger, the heat flow at well depth is typically not the same between adjacent time steps. To describe the change in heat flow over time, this application defines the heat flow increment corresponding to the nth time step as... ,in This represents the increase or decrease in heat flux per unit well depth in the current time step relative to the previous time step. If the heat transfer per unit well depth at a certain calculation node was 80 W / m in the previous time step, and it increases to 85 W / m in the current time step, then the heat flux increment at the corresponding node is 5 W / m.

[0025] Subsequently, based on the thermal response matrix established in step S31, the contribution of the heat flux increment at each historical time step to the current wellbore temperature is calculated. Specifically, each historical time step corresponds to a time-lag thermal response matrix. ,in This represents the time lag order from the current time step; a smaller time lag indicates that the time is closer to the current time, while a larger time lag indicates that the time is further away from the current time. For each historical time step, the heat response matrix with the corresponding time lag is first used. Corresponding historical heat flow increment Perform matrix multiplication to obtain the wellbore temperature rise generated by the historical heat flow at the current moment; then, accumulate the temperature rise results corresponding to all historical time steps in chronological order, and compare them with the initial wellbore temperature corresponding to the initial geothermal field. By superimposing the data, the wellbore temperature at the current time step can be obtained, and its expression is: As can be seen from the above recursive process, the wellbore temperature is essentially the result of the combined effect of the initial geothermal field and all historical heat flows. The more recent the historical heat flow is, the more obvious its temperature influence is usually. The more distant the historical heat flow is, although its influence gradually weakens, its temperature contribution will not disappear immediately due to the large heat capacity of the underground surrounding rock, but will continue to participate in the subsequent heat transfer process.

[0026] For example, after a deep geothermal heat exchanger has been operating continuously for a year, even if the operating conditions remain stable at a certain point, the heat injected underground over a long period has not yet fully dissipated, and there is still a certain amount of heat accumulation inside the surrounding rock. Therefore, the current wellbore temperature is still affected by the historical heat flow over several months or even longer, and will not return to the temperature state corresponding to the initial geothermal field. Using the aforementioned historical heat flow convolution recursive method can fully consider the heat accumulation and diffusion process in the underground surrounding rock, accurately describing the dynamic changes in wellbore temperature during the long-term operation of the deep geothermal heat exchanger. Compared to steady-state heat transfer calculation methods that only consider the current heat flow, this method is more consistent with the underground unsteady-state heat transfer mechanism and can effectively improve the accuracy of subsequent wellbore fluid temperature calculations and operating state predictions.

[0027] S33. Using the wellbore temperature as the boundary condition for the coupled heat transfer equation of the wellbore fluid, and combining it with the thermal resistance network between the wellbore fluid, the pipe wall and the surrounding rock, solve for the temperature distribution of the wellbore fluid along the pipe. Specifically, in step S33, after obtaining the wellbore temperature calculated in step S32, the coupled heat transfer relationship between the wellbore fluid and the surrounding rock is further established, and the temperature distribution along the wellbore fluid is solved. In this application, the wellbore temperature is used as the boundary condition for the wellbore fluid heat transfer calculation to describe the process of continuous heat transfer from the surrounding rock to the wellbore interior. Since the wellbore fluid exchanges heat with the wellbore wall and between different flow channels during its flow within the well, it is necessary to establish a coupled heat transfer model between the wellbore fluid, the wellbore wall, and the surrounding rock to accurately describe the dynamic changes in the temperature inside the wellbore along the depth and time directions.

[0028] This application uses a coaxial casing deep geothermal heat exchanger as an example. The wellbore fluid consists of two parts: annular fluid and central tube fluid. The circulating fluid first enters the annulus from the wellhead, flows along the wellbore towards the bottom of the well, and continuously absorbs heat released by the surrounding rock. After reaching the bottom of the well, the fluid reverses course, enters the central tube, and then flows back towards the wellhead, eventually returning to the surface. Therefore, during the entire heat exchange process, three heat transfer modes exist simultaneously inside the wellbore: radial heat conduction between the surrounding rock and the annular fluid, radial heat transfer between the annular fluid and the central tube fluid, and convective heat transfer generated by the fluid flowing axially along the wellbore. These various heat transfer processes are coupled with each other and jointly determine the temperature distribution of the fluid inside the wellbore.

[0029] To describe the heat transfer process described above, this application employs a thermal resistance network to establish the coupled heat transfer relationship within the wellbore. Specifically, a wellbore thermal resistance is established between the wellbore wall and the annular fluid. This represents the overall thermal resistance encountered during the transfer of heat from the surrounding rock through the wellbore to the annular fluid; a thermal resistance is set between the annular fluid and the fluid in the central tube. This represents the combined thermal resistance formed by heat conduction through the central pipe wall and convective heat transfer inside and outside the pipe. Therefore, the heat released from the surrounding rock first passes through the well wall thermal resistance. It is transferred to the annular fluid, and then through the thermal resistance. The heat is transferred to the fluid in the central tube, thus forming a continuous coupled heat transfer network between the surrounding rock, wellbore, annular fluid, and central tube fluid. The aforementioned thermal resistance parameters can be calculated based on the wellbore structural dimensions, material thermal conductivity, and fluid heat transfer coefficient, accurately reflecting the heat transfer capacity between different media within the wellbore. After establishing the thermal resistance network, coupled heat transfer equations for the annular fluid and central tube fluid are established according to the principle of fluid energy conservation. For the annular fluid, its temperature change is simultaneously affected by three factors: first, the heat carried by the fluid flowing in from the previous depth calculation node; second, the heat transferred through the wellbore via thermal resistance... The heat transferred to the annular fluid; and thirdly, the heat transferred through the annular fluid via thermal resistance. The heat transferred to the fluid in the central tube. For the fluid in the central tube, its temperature change is mainly affected by the heat flowing in from the fluid at the previous depth node and the heat transfer resistance of the annular fluid. The combined effect of heat transfer allows us to establish a coupled energy conservation relationship between different flow channels in the wellbore.

[0030] To improve the stability of numerical calculations, this application employs an implicit upwind discretization method to discretize the aforementioned heat transfer equations. Implicit discretization means that all unknown temperatures at the current time step participate in establishing the discretization equations, thus ensuring good numerical stability throughout the calculation process. Even with a large time step, the calculation results can be guaranteed to converge stably, minimizing numerical oscillations. Upwind discretization establishes spatial discretization relationships based on the actual fluid flow direction, meaning that the downstream node temperature is mainly influenced by the upstream node temperature. For example, when fluid flows along the annulus from the wellhead to the bottom, the annular fluid temperature at the j-th depth calculation point is mainly influenced by the upstream node temperature. The temperature of a single calculation point is determined by the heat carried by the fluid flowing into it. Similarly, when the fluid flows along the central pipe from the bottom of the well to the wellhead, its temperature is calculated node by node according to the flow direction of the central pipe. By adopting implicit upwind discretization, not only can the actual flow direction of the fluid be kept consistent with the numerical discretization direction, but the numerical oscillations that are prone to occur under convection-dominated conditions can also be effectively avoided, thus improving the accuracy of wellbore fluid temperature calculation.

[0031] After discretizing the heat transfer equations at each depth node, the discrete heat transfer equations for the annular fluid, the discrete heat transfer equations for the central tube fluid, and the wellhead inlet boundary conditions and bottom-hole return boundary conditions are combined to form a unified set of linear equations: ;in, This is the wellbore fluid coupling coefficient matrix, used to describe the coupling relationships between all unknown temperatures at the current time step. This coefficient matrix is ​​jointly determined by parameters such as the equivalent linear heat capacity of the annular and central pipe channels, the circulating fluid mass flow rate, the fluid specific heat capacity, the wellbore wall thermal resistance, the inter-fluid heat transfer resistance, the depth step size, and the time step size. Therefore, it can comprehensively reflect the heat transfer and fluid flow relationships between various computational nodes within the wellbore. Unknown quantity vector This column vector, representing all temperatures to be solved at the current time step, including all annular fluid temperatures from the wellhead to the bottom and all central tube fluid temperatures from the bottom to the wellhead, can be represented as: ; in, This represents the annular fluid temperature corresponding to the j-th depth calculation node at the n-th time step. This represents the fluid temperature in the center tube corresponding to the calculation node at the j-th depth in the n-th time step, and N represents the total number of calculation nodes after discretization in the well depth direction. The load vector on the right-hand side of the linear equation system. The load vector, used to represent the known heat input at each calculation node in the current time step, mainly comprises three parts: the first part is the heat transferred from the surrounding rock to the well fluid corresponding to the wellbore temperature calculated in step S32, reflecting the continuous heating effect of the underground surrounding rock on the well fluid; the second part is the known terms formed by boundary conditions such as the wellhead inlet temperature and the well bottom return connection conditions, used to ensure the continuity of fluid flow and temperature; the third part is the time term formed by the fluid temperature at each calculation node in the previous time step, used to reflect the fluid's own thermal inertia and the time accumulation effect. Therefore, the load vector comprehensively reflects the external heat input to the well fluid in the current time step and the thermal state in the previous time step, providing known conditions for solving the wellbore temperature field at the current moment. After obtaining the coefficient matrix... and load vector Subsequently, Gaussian elimination, LU decomposition, or iterative algorithms for large sparse matrices can be used to solve the aforementioned linear equations, obtaining the annular fluid temperature and central tube fluid temperature corresponding to all depth calculation nodes at the current time step, thus yielding the temperature distribution of the wellbore fluid along the well depth direction. Since this temperature distribution simultaneously considers factors such as surrounding rock heat transfer, fluid coupling heat transfer, axial convection of the fluid, and historical wellbore temperature variations, it can accurately reflect the actual temperature variation patterns inside the wellbore under different operating conditions of the deep geothermal heat exchanger. This provides a reliable data foundation for subsequent calculations of analytical heat transfer model outputs such as outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay indices.

[0032] S34. Determine the analytical heat transfer model output based on the temperature distribution along the wellbore fluid. The analytical heat transfer model output includes at least the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term heat decay index. Specifically, in step S34, after obtaining the temperature distribution of the wellbore fluid along the depth direction, the analytical heat transfer model output can be further determined. The analytical heat transfer model output is used to characterize the heat transfer performance and operating status of the deep geothermal heat exchanger under current operating conditions and serves as the physical prior information for subsequent physical-aided machine learning models. The analytical heat transfer model output includes at least the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term heat decay index.

[0033] The outlet temperature refers to the fluid temperature when the circulating fluid returns to the wellhead through the central pipe after completing downhole heat exchange. In this embodiment, the temperature of the calculated node corresponding to the wellhead central pipe in the wellbore fluid temperature distribution obtained by step S33 is the outlet temperature, i.e.: ; in, This represents the fluid temperature corresponding to the calculation node of the wellhead center pipe at the nth time step. The outlet temperature directly reflects the current heating or cooling capacity of the deep geothermal heat exchanger and is an important indicator for evaluating system performance. For example, during winter heating operation, a higher outlet temperature indicates that the circulating fluid absorbs more heat from the ground, and the system has a stronger heating capacity; however, with long-term operation, if the outlet temperature gradually decreases, it indicates that the heat from the surrounding rock is gradually decreasing, and the heat exchange capacity begins to decline.

[0034] The heat transfer per unit depth represents the heat transfer capacity per unit length of the wellbore along the depth direction, reflecting the contribution of different well sections to the overall heat transfer. Instantaneous heat transfer represents the total heat transfer power between the deep geothermal heat exchanger and the surrounding underground rock at the current operating moment. It can be calculated based on the inlet temperature, outlet temperature, mass flow rate, and specific heat capacity of the circulating fluid. When the system is in heating mode, the outlet temperature is higher than the inlet temperature, indicating that the circulating fluid absorbs heat from the surrounding rock; when the system is in cooling mode, it indicates that the circulating fluid releases heat to the underground. Instantaneous heat transfer directly reflects the system's heat transfer capacity under current operating conditions and is an important indicator for evaluating system operating efficiency. Long-term thermal decay describes the degree of decrease in heat transfer capacity of a deep geothermal heat exchanger after long-term continuous operation due to the continuous extraction or accumulation of heat from the surrounding rock. During long-term operation, the temperature of the surrounding rock near the wellbore gradually deviates from the initial geothermal field, and the local heat reservoir gradually changes, leading to a gradual decrease in outlet temperature and heat transfer. To quantitatively evaluate this long-term heat decay phenomenon, this embodiment can use the change ratio of the outlet temperature or instantaneous heat exchange after long-term operation relative to the initial operation stage as the long-term heat decay index.

[0035] In this embodiment, the aforementioned outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay index are all derived from the calculation results of the analytical heat transfer model established in steps S31 to S33, which can comprehensively reflect the current and long-term operating status of the deep geothermal heat exchanger. Specifically, the outlet temperature and instantaneous heat transfer reflect the current operating performance of the equipment, the heat transfer per unit well depth reflects the heat transfer distribution characteristics along the depth direction of the wellbore, and the long-term thermal decay index reflects the long-term evolution law of the surrounding rock thermal field. Using the output of the analytical heat transfer model as the physical prior input of the physical-aided machine learning model not only enables the machine learning model to fully utilize information on underground heat transfer mechanisms but also effectively compensates for the lack of field-measured samples, improving the accuracy, generalization ability, and long-term prediction stability of the model's predictions of operating status under different wellbore structures, formation conditions, and operating conditions. Thus, the machine learning model inherits the physical constraints of the analytical heat transfer model while possessing the nonlinear fitting capability of a data-driven model, thereby achieving complementary advantages between the physical model and the data model.

[0036] S4. Using the analytical heat transfer model, parameter calculations are performed under different wellbore structure parameters, operating condition parameters, and drilling exploration data to generate simulated operating data. Based on the simulated operating data, adaptive enhancement is performed to construct a simulated dataset. Specifically, in step S4, the analytical heat transfer model established in step S3 is used to perform parametric calculations under different combinations of wellbore structural parameters, drilling exploration data, and operating condition parameters. This yields a large amount of simulated operating data covering various engineering scenarios. Adaptive enhancement is then performed based on this simulated operating data, ultimately forming a simulated dataset for training a physical-aided machine learning model. Because on-site testing of deep geothermal heat exchangers is costly and has a long operating cycle, the amount of actual measured operating data obtainable is usually limited. Therefore, generating a large amount of physically constrained simulated operating data through the analytical heat transfer model can effectively compensate for the lack of actual measured samples, providing sufficient data support for subsequent machine learning model training.

[0037] In this embodiment, parametric calculation refers to automatically calling the analytical heat transfer model to repeatedly calculate the operating results corresponding to each working condition according to a pre-set combination of multiple input parameters, so that the simulation data covers the largest possible parameter space. The input parameters include wellbore structural parameters, drilling exploration data, and operating condition parameters. Wellbore structural parameters include, for example, well depth, well diameter, well wall radius, length of the wellbore heat exchange section, inner diameter of the central tube, outer diameter of the central tube, annular dimensions, and casing dimensions; drilling exploration data includes the top and bottom depths of each formation, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradient; operating condition parameters include inlet temperature, circulation flow rate, heat load, and operating time. For example, in a deep geothermal heat exchanger project, the well depth can be set to vary from 2500m to 4000m, the inlet temperature from 5℃ to 30℃, the circulation flow rate from 10kg / s to 40kg / s, the formation thermal conductivity from 1.8W / (m·K) to 4.5W / (m·K), and the geothermal gradient from 20℃ / km to 40℃ / km. The analytical heat transfer model calculates the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay index under different parameter combinations, thus generating simulated operating data.

[0038] Specifically, S4 in this embodiment includes: S41. Based on the wellbore structure parameters, operating condition parameters and drilling exploration data, determine the input parameters and the value range of each input parameter for the analytical heat transfer model. The initial simulation conditions are generated by combining Latin hypercube sampling, Sobol low-difference sequence sampling and equally spaced layered sampling. Latin hypercube sampling is used to ensure that each input parameter is evenly distributed across the entire value range, so that every parameter interval can be sampled, avoiding sample clustering that can occur with random sampling. For example, when the well depth range is 2500m to 4000m, it can be divided into 100 intervals, and one sample can be randomly drawn from each interval to cover the entire well depth range.

[0039] Sobol low-discrepancy sequence sampling is used for uniform filling of the multidimensional parameter space. Since the input parameters of deep geothermal heat exchangers typically include dozens of variables, ordinary random sampling alone can easily generate a large number of blank areas. Sobol sequences, however, can make the sample distribution in the high-dimensional parameter space more uniform. For example, for the five-dimensional parameter space of well depth, inlet temperature, circulation flow rate, thermal conductivity, and geothermal gradient, Sobol sequences can ensure that the samples uniformly cover the entire five-dimensional hypercube, improving the training quality of subsequent surrogate models. Equal-interval stratified sampling is mainly used to ensure that there are sufficient samples near the boundaries of important parameters. For example, the inlet temperature is stratified at 5°C intervals, and the circulation flow rate is stratified at 5 kg / s intervals. Each layer generates a certain number of operating conditions, ensuring that different operating areas have representative samples. In this embodiment, the above three sampling methods are combined: Latin hypercube sampling ensures uniform marginal distribution of each variable, Sobol sampling improves the coverage of the high-dimensional space, and equal-interval stratified sampling ensures that important operating conditions are covered, thereby generating initial simulation conditions covering the entire parameter space. For example, 5,000 initial simulation conditions can be generated in the end, and each condition corresponds to a unique set of wellbore structural parameters, drilling exploration data and operating condition parameters.

[0040] S42. Input the initial simulation conditions into the analytical heat transfer model to obtain the analytical heat transfer model output corresponding to each initial simulation condition, and generate the corresponding initial simulation operation data. In step S42, the initial simulation conditions obtained in step S41 are input into the analytical heat transfer model established in step S3 for calculation, resulting in the corresponding analytical heat transfer model output and generating initial simulation operation data. In this embodiment, each set of simulation operation data includes model input parameters and model output parameters. The input part includes wellbore structure parameters, drilling exploration data, and operating condition parameters; the output part includes outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term heat decay index.

[0041] For example, for one set of simulation conditions: well depth 3500m; geothermal gradient 32℃ / km; circulation flow rate 25kg / s; inlet temperature 12℃; thermal conductivity 3.2W / (m·K). The analytical heat transfer model calculates: outlet temperature 38.4℃; instantaneous heat transfer 2.75MW; heat transfer per unit well depth 785W / m; long-term thermal decay index 0.082. These inputs and outputs together constitute a single set of simulation operation data. After completing the calculations for all initial simulation conditions sequentially, the initial simulation operation dataset can be formed.

[0042] S43. Based on the output of the analytical heat transfer model, determine the physical response sensitivity of each input parameter of the analytical heat transfer model, and construct a sampling scoring function by combining the prediction uncertainty of the surrogate model and the boundary distance of the operating condition. Based on the sampling scoring function, perform adaptive sampling on the parameter space to obtain the new simulated operating condition. In one embodiment of this application, the sampling score value corresponding to each candidate operating condition in the parameter space is calculated, and each candidate operating condition is sorted according to the sampling score value. The candidate operating conditions with a sampling score value higher than a preset threshold or the top preset number of candidate operating conditions are selected as new simulation conditions to achieve adaptive sampling of the parameter space.

[0043] The sampling scoring function is expressed as follows: ; in, This indicates the uncertainty in the surrogate model prediction. This represents the physical response sensitivity corresponding to the output of the analytical heat transfer model. Indicates the distance to the working condition boundary. Indicates indicators under extreme operating conditions. These represent the corresponding weights; For sampling scoring functions; In detail, in step S43, in order to further improve the coverage of complex working conditions and boundary working conditions of the simulation dataset, this embodiment does not continue to randomly generate samples, but rather performs adaptive enhancement based on existing simulation running data.

[0044] First, the physical response sensitivity corresponding to each input parameter is determined based on the output of the analytical heat transfer model. Physical response sensitivity refers to the degree to which the output of the analytical heat transfer model responds to changes in each input parameter, and can be calculated using local perturbation analysis. For example, keeping other parameters constant, increasing the well depth by 1% and recalculating the analytical heat transfer model output, the corresponding sensitivity is calculated based on changes in outlet temperature, instantaneous heat transfer, etc. If the outlet temperature changes significantly, it indicates a high physical response sensitivity corresponding to the well depth; if the change is small, it indicates a low sensitivity for that parameter.

[0045] Subsequently, a surrogate model is established. The surrogate model can be a Gaussian process regression model, a random forest model, or a neural network model, used to approximate the mapping relationship between the input and output of the analytical heat transfer model. Since the surrogate model is computationally much faster than the analytical heat transfer model, it can quickly estimate the prediction results and prediction uncertainties at various locations throughout the entire parameter space.

[0046] Among them, the uncertainty of the surrogate model prediction This is used to indicate the reliability of the surrogate model's prediction results for the current parameter point. Generally, densely sampled regions have lower prediction uncertainty, while sparsely sampled regions have higher prediction uncertainty. Therefore, prioritizing the generation of samples in regions with higher prediction uncertainty can improve the overall coverage of the parameter space.

[0047] Operating boundary distance This indicates the distance of the current parameter point from the parameter value boundary. For example, when the allowable well depth range is 2500m to 4000m, the point near 3900m is closer to the boundary, with a smaller boundary distance; the point near 3000m is farther from the boundary, with a larger boundary distance. Since extreme conditions in engineering design usually occur near the boundary, the sample size needs to be increased near the boundary.

[0048] Extreme operating condition indicators This is used to describe whether the current parameters are close to extreme operating conditions such as high temperature, high flow rate, and large well depth. For example, when the local temperature gradient exceeds 38℃ / km, the well depth exceeds 3900m, or the inlet temperature is below 5℃, the current operating condition can be determined to be in an extreme operating zone, and a higher extreme operating condition index can be assigned.

[0049] Subsequently, a sampling scoring function is constructed based on the above indicators: ; in, These represent the corresponding weights, which can be determined based on experience or cross-validation. For example, they can be set as follows: =0.35, =0.30, =0.20, =0.15; For each candidate operating condition in the parameter space, the corresponding score is calculated. The higher the score, the more worthwhile the operating condition is to continue the analytical heat transfer model calculation. For example: Candidate operating condition A: prediction uncertainty 0.82; sensitivity 0.75; boundary distance 0.60; extreme condition index 0.20; then: S = 0.35 × 0.82 + 0.30 × 0.75 + 0.20 × 0.60 + 0.15 × 0.20 = 0.662; Candidate operating condition B: prediction uncertainty 0.45; sensitivity 0.35; boundary distance 0.40; extreme condition index 0.10; then: S = 0.338; Since the score of operating condition A is significantly higher than that of operating condition B, operating condition A is selected first for analytical heat transfer model calculation.

[0050] Therefore, the adaptive sampling described in this embodiment is not random sampling, but rather automatically selects the most valuable working conditions from the parameter space for supplementary calculation based on the sampling scoring function. This ensures that the newly added simulated working conditions are mainly concentrated in the uncertain region of model prediction, the region of drastic change in physical response, and the region of extreme engineering working conditions. This reduces the number of calculations of the analytical heat transfer model while ensuring the quality of the simulation data and improves the efficiency of dataset construction.

[0051] S44. Input the newly added simulated operating conditions and the preset extreme operating conditions into the analytical heat transfer model to generate corresponding newly added simulated operating data, and combine the newly added simulated operating data with the initial simulated operating data to form the simulation dataset.

[0052] Specifically, in step S44, the newly added simulated operating conditions obtained in step S43 and the pre-set extreme operating conditions are input into the analytical heat transfer model for calculation to obtain the corresponding newly added simulated operating data. Among them, the pre-set extreme operating conditions may include engineering extreme conditions such as maximum well depth, maximum circulation flow rate, minimum inlet temperature, maximum geothermal gradient, and minimum thermal conductivity, which are used to further enhance the machine learning model's ability to learn extreme operating states.

[0053] For example, based on an initial set of 5000 simulated operating conditions, 1500 high-value simulated operating conditions are added using a sampling scoring function, followed by 500 extreme operating conditions, ultimately resulting in 7000 sets of simulated operating data. Finally, the newly added simulated operating data is combined with the initial simulated operating data to construct the final simulation dataset. This simulation dataset not only covers normal operating conditions but also regions with drastic changes in physical response, parameter space boundary regions, and extreme operating regions. This allows subsequent physics-assisted machine learning models to learn more comprehensive input-output mapping relationships, improving model prediction accuracy and long-term operational stability.

[0054] S5. Use the simulated dataset and the measured running data to train the physics-assisted machine learning model and obtain the trained physics-assisted machine learning model. The physics-assisted machine learning model includes a static encoder, a temporal encoder, a physical prior encoder, and multiple decoders, wherein: The static encoder is used to encode wellbore structural parameters and drilling exploration data; The timing encoder is used to encode timing data corresponding to inlet temperature, flow rate, heat load, start / stop status, and historical operating status. The physical prior encoder is used to encode the output of the parsed heat transfer model and the response kernel features extracted from the thermal response matrix; The features output by the static encoder, the temporal encoder, and the physical prior encoder are fused to form a unified latent variable: ;in, Represents a unified latent variable. This represents the static coding features output by the static encoder. This represents the timing coding characteristics output by the timing encoder. F represents the physical prior encoded features output by the physical prior encoder, and F() represents the feature fusion function; The unified latent variable Multiple decoders are input, and the outlet temperature, heat exchange capacity, long-term heat decay index and prediction deviation correction amount are output respectively to obtain the corrected deep geothermal heat exchanger operation prediction model. The multiple decoders include an outlet temperature decoder, a heat exchanger decoder, a long-term thermal decay decoder, and a deviation correction decoder. The outlet temperature decoder is used to output the outlet temperature of the deep geothermal heat exchanger. The heat exchanger decoder is used to output the heat exchanger of the deep geothermal heat exchanger. The long-term thermal decay decoder is used to output the long-term thermal decay index of the deep geothermal heat exchanger. The deviation correction decoder is used to learn and analyze the prediction deviation between the output of the heat transfer model and the measured operation results, and output the corresponding prediction deviation correction amount. The predicted deviation correction amount is fused with the output of the analytical heat transfer model to obtain the corrected prediction result for the deep geothermal heat exchanger operation, which is expressed as follows: ; in, This represents the output of the analytical heat transfer model. This indicates a deviation correction decoder. The corrected operating prediction results for the deep geothermal heat exchanger include outlet temperature, heat exchange capacity, and long-term heat decay index.

[0055] The training objective function of the physics-assisted machine learning model is expressed as: ; in, This represents the prediction error corresponding to the simulated dataset; This represents the prediction error corresponding to the measured operating data; This represents the time series smoothing constraint error of the time series data corresponding to the operating condition parameters; This represents the domain deviation correction error between the simulated dataset and the measured operating data; , , and These represent the weighting coefficients of the corresponding loss terms; For training the objective function; When the training objective function converges to a preset loss threshold, or when the change in the training objective function within a preset number of consecutive training rounds is less than a preset convergence threshold, the training of the physical-assisted machine learning model is determined to be complete, and a trained physical-assisted machine learning model is obtained.

[0056] Next, we will explain step S5 in detail with specific examples: In step S5, the simulated dataset obtained in step S4 and the measured operational data collected on-site are used together to train a physical-aided machine learning model, resulting in a trained prediction model for the operating status of the deep geothermal heat exchanger. The simulated dataset originates from large-scale parameterized calculations of the analytical heat transfer model, which can fully cover different wellbore structural parameters, different drilling exploration data, and different operating conditions, exhibiting good physical consistency. The measured operational data, on the other hand, comes from actual engineering operations, reflecting complex disturbances, equipment errors, and environmental changes in real-world engineering projects. Therefore, this application uses both simulated and measured operational data for model training, enabling the machine learning model to learn both the underground heat transfer mechanism and the actual operating patterns in real-world engineering projects, thus balancing physical reliability and prediction accuracy. The measured operational data includes operating condition parameters and corresponding actual operating results; wherein the operating condition parameters include at least inlet temperature, circulation flow rate, heat load, and operating time, and the actual operating results include at least one or more of the following: outlet temperature, heat exchange rate, and long-term heat decay index.

[0057] For example, for a deep geothermal heat exchanger project, step S4 can generate approximately 8,000 sets of simulated operation data, while simultaneously collecting approximately 500 sets of field operation data. The simulated operation data covers different operating conditions, such as well depths of 2500m to 4000m, geothermal gradients of 20℃ / km to 40℃ / km, and circulation flow rates of 10kg / s to 40kg / s. The field data reflects the actual outlet temperature, heat exchange, and heat decay during long-term operation of the equipment. Both types of data participate in model training, effectively improving the model's generalization ability under unknown operating conditions. In this embodiment, the physical-aided machine learning model includes a static encoder, a temporal encoder, a physical prior encoder, and multiple decoders. Each encoder extracts different types of data features, which are then fused together, and finally, multiple decoders predict different operating indicators. The static encoder is used to extract static features corresponding to wellbore structural parameters and drilling exploration data. Static data refers to data that remains unchanged or changes very slowly during a single prediction process. For example, wellbore structural parameters such as well depth, well diameter, well wall radius, inner diameter of the central tube, outer diameter of the central tube, annulus dimensions, casing dimensions, and length of the heat exchange section, as well as drilling exploration data such as the top and bottom depths of each formation, layer thickness, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradient, do not change during equipment operation and are therefore considered static information. For instance, for a deep geothermal heat exchanger, the following static input can be established: well depth 3500m; well wall radius 0.108m; outer diameter of the central tube 0.114m; inner diameter of the central tube 0.098m; 8 formations; thermal conductivity of the first layer 2.4W / (m·K); thermal conductivity of the second layer 3.1W / (m·K); ...; corresponding geothermal gradients for each layer are 24℃ / km, 28℃ / km, and 32℃ / km, respectively. All the above parameters are arranged in a uniform order to form a static input vector, which is then input into the static encoder. Static encoders can employ multi-layer fully connected neural networks for encoding, transforming the original high-dimensional static parameters into fixed-dimensional static encoded features through layer-by-layer nonlinear mapping. For example, the original static parameters, which consist of 40 dimensions, can be mapped to 64-dimensional static features after being processed by a static encoder, thereby preserving information on the influence of wellbore structure and formation conditions on the heat transfer process.

[0058] A time-series encoder is used to extract the dynamic characteristics of operating conditions that change over time. Unlike wellbore structural parameters, inlet temperature, circulation flow rate, heat load, start-up and shutdown status, and historical operating status all change continuously over time, thus belonging to time-series data. For example, if a heat exchanger operates continuously for 24 hours, recording operating parameters once per hour, the inlet temperature can be represented as: 12℃, 12.5℃, 13℃, 13.4℃, ...; the circulation flow rate is: 20kg / s, 22kg / s, 21kg / s, 23kg / s, ...; and the historical outlet temperature is: 35℃, 35.4℃, 35.8℃, 36.1℃, ...; the above data together constitute the time-series input. In this embodiment, the time-series encoder can be implemented using a Long Short-Term Memory (LSTM) network, a Gated Recurrent Unit (GRU), or a Transformer network. Because deep geothermal heat exchangers have significant thermal inertia, the current operating status is usually influenced by the operating status several hours or even days ago. Therefore, using a network with long-term time-series memory capabilities can effectively learn the dynamic changes during operation. Timing encoder output ;in, This represents time-series coding features. For example, 48 consecutive hours of operational data can be encoded into 128-dimensional dynamic features to describe the heat accumulation patterns formed by changes in operational status over time.

[0059] The physical prior encoder is used to extract the physical mechanism features corresponding to the output of the analytical heat transfer model and the thermal response matrix. Unlike ordinary data-driven models, this application does not directly use operational data for prediction. Instead, it first uses the analytical heat transfer model to obtain calculation results with clear physical meaning, such as outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay index. At the same time, it further uses the thermal response matrix established in step S31 to extract the underground heat transfer response features.

[0060] The analytical heat transfer model output reflects the underground heat conduction pattern, while the thermal response matrix reflects the thermal coupling relationship between different underground locations. For example, for 150 well depth calculation nodes and 300 time steps, a thermal response matrix with a size of 150×150×300 can be obtained.

[0061] Since the thermal response matrix has a high dimension, it is not suitable to directly input it into a neural network. Therefore, feature extraction is performed first. For example, average pooling, convolution dimensionality reduction, or principal component analysis are used to extract response kernel features, compressing the high-dimensional thermal response matrix into a fixed-length physical response vector.

[0062] For example, a 64-dimensional response kernel feature can be extracted, which reflects the thermal inertia, heat propagation velocity, and long-term heat accumulation characteristics of the underground surrounding rock. Subsequently, the output of the analytical heat transfer model and the response kernel feature are input together into the physical prior encoder to obtain the physical prior encoded features. Since this feature originates from the heat transfer mechanism, it can provide physical constraints for machine learning models, preventing the models from producing predictions that violate the laws of heat transfer.

[0063] After the three encoders have been computed, feature fusion is required. In this embodiment, statically encoded features are fused. Time-series coding features and physical prior coding features By merging them, a unified latent variable can be formed: F() represents the feature fusion function, which can be implemented by feature concatenation, weighted fusion, or attention fusion.

[0064] Since deep geothermal heat exchangers need to predict multiple operating indicators simultaneously, this application adopts a multi-task learning approach to establish multiple decoders, including an outlet temperature decoder, a heat exchange rate decoder, a long-term heat decay decoder, and a deviation correction decoder.

[0065] Among them, the outlet temperature decoder is used to predict the wellhead outlet temperature; the heat transfer decoder is used to predict the current heat transfer capacity of the system; the long-term heat decay decoder is used to predict the performance change trend after long-term operation; and the deviation correction decoder is used to learn and analyze the systematic errors between the heat transfer model and the actual engineering.

[0066] S6. Obtain the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted. Based on the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted, obtain the output of the analytical heat transfer model corresponding to the operating condition to be predicted. Input the wellbore structural parameters, drilling exploration data, operating condition parameters, and the output of the analytical heat transfer model into the trained physical-aided machine learning model to obtain the predicted operating status of the deep geothermal heat exchanger.

[0067] Specifically, when predicting the operating status of a deep geothermal heat exchanger under a specific condition, the wellbore structural parameters, drilling exploration data, and operating parameters for that condition are first obtained. The wellbore structural parameters describe the heat exchanger's geometry, including well depth, well diameter, well wall radius, wellbore heat exchange section length, central tube inner diameter, central tube outer diameter, annular inner diameter, annular outer diameter, casing inner diameter, casing outer diameter, central tube wall thickness, and casing wall thickness. These parameters determine the fluid flow space, heat exchange area, and thermal resistance network structure within the wellbore, and are crucial fundamental parameters affecting heat exchange capacity. For example, a coaxial deep geothermal heat exchanger with a well depth of 2500m, a well diameter of 216mm, a central tube outer diameter of 114mm, and an annular outer diameter of 177mm has a significantly larger heat exchange area than a heat exchanger with a well depth of 1800m, thus its long-term heat exchange capacity is typically improved.

[0068] After obtaining drilling exploration data, it is used to reflect the true geological conditions of the surrounding underground rock. Drilling exploration data includes at least the top and bottom depths of each stratum, layer thickness, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradient. For example, in a certain project, the 0–600m layer is sandstone, with a thermal conductivity of approximately 2.4 W / (m·K) and a geothermal gradient of approximately 2.6℃ / 100m; the 600–1700m layer is limestone, with a thermal conductivity of approximately 3.3 W / (m·K) and a geothermal gradient of approximately 3.1℃ / 100m; and the 1700–2500m layer is granite, with a thermal conductivity of approximately 3.8 W / (m·K) and a geothermal gradient of approximately 3.5℃ / 100m. Because the thermal properties of different strata vary significantly, their heat transfer capabilities also differ considerably. The aforementioned drilling exploration data can accurately reflect the characteristics of heterogeneous underground strata, providing reliable geological input conditions for analyzing heat transfer models.

[0069] Simultaneously, the operating parameters corresponding to the predicted operating conditions are also acquired. These operating parameters describe the current actual operating status of the equipment, such as inlet temperature, circulation flow rate, heat load, operating time, and start / stop status. For example, under winter heating conditions, the inlet temperature may be 18℃, the circulation flow rate may be 25kg / s, the heat load may be 420kW, and the continuous operating time may be 180 days; while under summer cooling conditions, the inlet temperature may be adjusted to 30℃, the circulation flow rate may be adjusted to 20kg / s, and the corresponding heat exchange process will also change. Therefore, different operating conditions correspond to different wellbore heat transfer processes.

[0070] After obtaining the above parameters, the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted are input into the analytical heat transfer model established in steps S1 to S3. The analytical heat transfer model first constructs an initial geothermal field based on the drilling exploration data, then calculates the wellbore temperature by combining the layered thermal response relationship, thermal response matrix, and historical heat flux convolution recursively, and then solves the temperature distribution along the wellbore fluid based on the wellbore temperature, finally obtaining the analytical heat transfer model output corresponding to the current operating condition.

[0071] The analytical heat transfer model outputs physical prior information that conforms to the constraints of the underground heat transfer mechanism. This information includes at least the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term heat decay index. For example, for the aforementioned 2500m deep geothermal heat exchanger, under a certain predicted operating condition, the analytical heat transfer model calculates an outlet temperature of 36.8℃, a heat transfer per unit well depth of 108W / m, an instantaneous heat transfer of 270kW, and a long-term heat decay index of 0.87, indicating that after long-term operation, its heat transfer capacity will retain approximately 87% of its initial operating capacity. Since these results are derived from the analytical heat transfer model, the prediction results can be guaranteed to meet the basic heat transfer laws.

[0072] Subsequently, wellbore structural parameters, drilling exploration data, operating condition parameters, and the output of the analytical heat transfer model are all input into the pre-trained physical-aided machine learning model. Specifically, the wellbore structural parameters and drilling exploration data are input into a static encoder to extract static features reflecting the wellbore structure and geological conditions; the operating condition parameters are input into a time-series encoder to extract time-series features corresponding to changes in operating status; and the output of the analytical heat transfer model is input into a physical prior encoder to extract physical prior features reflecting the underground heat transfer patterns. After each of the three encoders completes its information extraction, a unified latent variable is formed through a feature fusion function. The unified latent variable synthesis includes wellbore structure information, formation thermal property information, operating status information, and physical model calculation results, thus providing a comprehensive description of the current operating status of deep geothermal heat exchangers.

[0073] Next, the hidden variables will be unified. Multiple decoders are used for joint prediction. The outlet temperature decoder outputs the predicted outlet temperature; the heat transfer decoder outputs the predicted heat transfer; the long-term heat decay decoder outputs the long-term heat decay index; and the deviation correction decoder outputs the deviation correction amount corresponding to the prediction error of the analytical heat transfer model. Since the analytical heat transfer model inevitably contains certain simplification errors, such as not fully considering complex factors like underground fissure seepage, groundwater convection heat transfer, or construction errors, the deviation correction decoder uses a large amount of measured operational data to learn the differences between the physical model and the actual operating results. For example, if the analytical heat transfer model predicts an outlet temperature of 36.8℃, and the deviation correction decoder outputs a correction amount of +0.9℃, the final corrected predicted outlet temperature is 37.7℃; if the analytical heat transfer model predicts an instantaneous heat transfer of 270kW, and the deviation correction amount is +8kW, the final predicted heat transfer is 278kW. Through deviation correction, the prediction accuracy can be further improved, making the prediction results closer to the actual operating conditions.

[0074] Finally, the physics-assisted machine learning model outputs predictions of the operating status of the deep geothermal heat exchanger. These predictions include at least the outlet temperature, heat exchange capacity, and long-term heat decay index. They can also be expanded to include wellbore temperature variation trends, fluid temperature distribution along the wellbore, heat exchange capacity distribution per unit well depth, and performance variation trends over future operating cycles, depending on specific application requirements. For example, it can predict that the outlet temperature will decrease from 37.7℃ to 36.9℃ within the next 30 days, and the long-term heat decay index will decrease from 0.87 to 0.85. This allows for advance assessment of heat exchanger performance trends, providing a basis for decision-making regarding operation scheduling, heat load allocation, and equipment maintenance.

[0075] Therefore, this application does not directly rely on machine learning models for prediction during the prediction stage. Instead, it first uses an analytical heat transfer model to calculate physical prior information that conforms to the laws of underground heat transfer. Then, it combines wellbore structural parameters, drilling exploration data, and operating condition parameters with the trained physical-aided machine learning model, achieving a deep integration of physical mechanisms and data-driven prediction. Compared with traditional analytical model prediction methods, it can effectively correct prediction errors caused by model simplification. Compared with pure data-driven prediction methods, it can achieve high prediction accuracy without relying on a large number of measured samples. It still has good generalization ability, prediction stability, and engineering application value under complex strata, small sample conditions, and long-term operational prediction conditions.

[0076] In some embodiments of this application, the method further includes: S7, acquiring newly measured operating data generated during the operation of the deep geothermal heat exchanger, and using the trained physical-aided machine learning model to predict the operating conditions corresponding to the newly measured operating data; when the prediction error between the prediction result and the measured result corresponding to the newly measured operating data exceeds a preset threshold, generating new simulated operating data using the operating conditions corresponding to the newly measured operating data, and incrementally updating the physical-aided machine learning model using the new simulated operating data and the new measured operating data.

[0077] Specifically, S7 includes: S71. Acquire newly measured operating data generated during the operation of the deep geothermal heat exchanger, and use the trained physical-aided machine learning model to predict the operating conditions corresponding to the newly measured operating data; compare the prediction results with the corresponding measured results, and calculate the prediction error; specifically, step S71 first acquires newly measured operating data generated during the operation of the deep geothermal heat exchanger. The newly measured operating data refers to the new operating data continuously collected during actual operation after the model training is completed, which can be acquired in real time by the on-site monitoring system. For example, operating data can be collected every 10 minutes, 30 minutes, or 1 hour, or statistically collected daily or weekly. In this embodiment, the newly measured operating data includes at least the inlet temperature, outlet temperature, circulation flow rate, heat load, operating time, wellhead pressure, bottom hole pressure, well wall temperature (if monitoring conditions are available), and heat exchange under the corresponding operating conditions.

[0078] For example, after two months of continuous operation, the on-site monitoring system acquired a new set of measured operating data for a deep geothermal heat exchanger. This data included an inlet temperature of 17.8℃, a circulation flow rate of 24.5 kg / s, a heat load of 390 kW, a continuous operating time of 1260 h, a measured outlet temperature of 36.9℃, and a measured instantaneous heat exchange of 275 kW. Subsequently, the wellbore structural parameters, drilling exploration data, operating condition parameters, and the output of the analytical heat transfer model corresponding to these new operating conditions were input into a pre-trained physical-aided machine learning model for prediction, yielding the corresponding prediction results. For example, the model predicted an outlet temperature of 35.8℃ and a predicted heat exchange of 265 kW.

[0079] Subsequently, the model prediction results are compared with the measured results in the corresponding newly added measured operating data, and the prediction error is calculated. The prediction error can be calculated using methods such as absolute error, relative error, mean square error (MSE), root mean square error (RMSE), or mean absolute percentage error (MAPE). For example, in this embodiment, the mean relative error is used as the evaluation index, so the outlet temperature prediction error is approximately 3.0%, and the heat exchange prediction error is approximately 3.6%. The above-mentioned multiple prediction indicators can also be comprehensively evaluated according to preset weights to obtain the comprehensive prediction error under the current operating conditions, which is used to determine whether the model still meets the prediction requirements.

[0080] S72. When the prediction error exceeds a preset threshold, a local disturbance operating condition is constructed centered on the operating condition corresponding to the newly added measured operating data. The local disturbance operating condition is input into the analytical heat transfer model established in step S1, and calculations are performed according to the thermal response matrix calculation, historical heat flux convolution recursion, and wellbore fluid coupling heat transfer solution process in step S3 to generate corresponding new simulated operating data. In step S72, when the prediction error exceeds the preset threshold, it indicates that the current operating condition has exceeded the range covered by the original training data of the model, or the underground heat transfer environment has changed significantly. At this time, the model incremental update process is initiated. The preset threshold can be set according to engineering requirements. For example, when the outlet temperature prediction error exceeds 2%, the heat exchange prediction error exceeds 5%, the long-term heat decay index prediction error exceeds 3%, or the comprehensive prediction error exceeds 4%, it is considered that a model update is required.

[0081] Once the update conditions are met, the operating condition corresponding to the newly added measured operating data is used as the center, and local disturbance operating conditions are constructed in its vicinity. Local disturbance operating conditions refer to a series of new operating conditions formed by small-scale disturbances to one or more operating parameters near the newly added measured operating condition, enabling the model to learn a richer data distribution in the vicinity of the current operating condition. For example, when the newly added operating condition corresponds to an inlet temperature of 18℃, a circulation flow rate of 25kg / s, and a heat load of 400kW, disturbance ranges such as inlet temperature ±1℃, circulation flow rate ±2kg / s, and heat load ±30kW can be set to form several new operating conditions. For example, different combinations such as inlet temperatures of 17℃, 18℃, and 19℃, circulation flow rates of 23kg / s, 25kg / s, and 27kg / s, and heat loads of 370kW, 400kW, and 430kW can be generated, thus forming dozens or even hundreds of local disturbance operating conditions.

[0082] Subsequently, the aforementioned local disturbance operating conditions are input into the analytical heat transfer model established in step S1, and the corresponding calculations are completed according to the calculation process in step S3. Specifically, firstly, the initial geothermal field is re-established based on the operating parameters after the disturbance; then, the thermal response matrix is ​​recalculated using the layered thermal response relationship established in step S2; next, the wellbore temperature is calculated recursively based on the historical heat flux convolution; and further, the coupled heat transfer equation of the wellbore fluid is solved to obtain the temperature distribution along the wellbore fluid. Finally, the analytical heat transfer model outputs, such as the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay index under the corresponding operating conditions, are calculated. The above calculation results, together with the corresponding input parameters, constitute the newly added simulated operating data.

[0083] For example, for the aforementioned local disturbance conditions, 100 new sets of simulated operating data can be generated, including outlet temperatures, heat exchange rates, and long-term heat decay indices corresponding to different inlet temperatures, different circulation flow rates, and different heat loads, thereby enriching the data samples near the current operating conditions.

[0084] S73. The newly added simulated running data and the newly added measured running data are used together as incremental training data to incrementally update the physical-assisted machine learning model. Among them, multiple decoders or the output layers corresponding to multiple decoders are updated first to correct the model prediction bias. In step S73, the newly added simulated running data and the newly added measured running data are combined to form incremental training data, which is used to incrementally update the physical-aided machine learning model. Since the model's front-end encoder has already learned relatively stable wellbore structural characteristics, formation thermal properties, and heat transfer mechanisms, this application prioritizes updating the output layer parameters of multiple decoders or multiple decoders, while keeping most of the parameters of the static encoder, temporal encoder, and physical prior encoder unchanged. This can quickly correct model prediction biases while avoiding the increased computational load and historical knowledge forgetting problems caused by retraining the entire network.

[0085] For example, when the model detects that the current outlet temperature prediction is consistently lower than expected by about 0.8℃, the deviation correction decoder can quickly learn the new deviation pattern and gradually increase the corresponding deviation correction amount from 0.2℃ to 1.0℃, thereby automatically correcting the analytical heat transfer model output and making the final prediction result closer to the actual measurement result. Similarly, for the problem of consistently low heat transfer prediction, the heat transfer decoder can also quickly complete parameter correction with a small number of new samples, improving the model prediction accuracy.

[0086] S74. When the cumulative number of newly added measured operating data reaches a preset number, the newly added simulated operating data and the newly added measured operating data are used together to fine-tune the physical-assisted machine learning model, complete the incremental update of the physical-assisted machine learning model, and use the incrementally updated physical-assisted machine learning model to predict the operating status of subsequent working conditions.

[0087] In step S74, to avoid parameter oscillations caused by frequent model updates, this application further sets up a cumulative update mechanism. When the cumulative number of newly added measured operational data reaches a preset quantity, the newly obtained simulated operational data and the newly added measured operational data are used together to jointly fine-tune the entire physical-aided machine learning model. The preset quantity can be set according to engineering needs, for example, joint fine-tuning can be initiated after accumulating 200, 500, or 1000 sets of newly added measured operational data. During the joint fine-tuning process, not only are the parameters of multiple decoders updated, but some network parameters of the temporal encoder and physical prior encoder can also be appropriately adjusted, allowing the model to gradually adapt to new underground heat transfer patterns and operating environments without retraining using all historical data, thus significantly reducing training costs. For example, after a deep geothermal heat exchanger has been running continuously for one year, a total of 500 sets of newly added measured operational data and 1000 sets of corresponding newly generated simulated operational data have been obtained, and the system automatically initiates joint fine-tuning. After fine-tuning, the average error in outlet temperature prediction decreased from 3.2% to 0.9%, the average error in heat exchange prediction decreased from 4.8% to 1.5%, and the error in long-term heat decay index prediction decreased from 5.1% to 1.8%. Subsequently, the physical-aided machine learning model with incremental updates replaced the original model and continued to be used for predicting the operating status of subsequent operating conditions. This enabled the model to continuously learn new operating rules as the equipment operates over a long period of time, maintaining high prediction accuracy, good generalization ability, and strong adaptability to complex operating conditions.

[0088] In summary, this application establishes an online incremental update mechanism driven by newly added measured operating data. When the model prediction error exceeds a preset threshold, it constructs a local disturbance operating condition using the current operating conditions and generates new simulated operating data in conjunction with the analytical heat transfer model. Then, it incrementally updates the physical-aided machine learning model in conjunction with the newly added measured operating data. This enables the model to continuously learn new characteristics, new patterns, and new operating conditions generated during the long-term operation of the deep geothermal heat exchanger. This not only avoids the decline in prediction performance caused by changes in the operating environment, but also fully leverages the advantages of combining the analytical heat transfer model and the machine learning model. While ensuring physical consistency, it improves the long-term prediction accuracy, stability, and reliability of the model in engineering applications.

[0089] Preferably, in some embodiments of this application, S31 specifically includes: S311. Based on the layered thermal response relationship, the multiple sub-paths corresponding to the finite linear heat source integral path are respectively taken as finite linear heat source segments in the corresponding strata, and the transient thermal response under the corresponding time lag is calculated using the equivalent thermal conductivity and equivalent thermal diffusivity corresponding to each sub-path. Specifically, step S311 first uses the layered thermal response relationship established in step S2 to divide the finite linear heat source integral path into multiple sub-paths, which are then used as finite linear heat source segments in the corresponding formations. Here, a finite linear heat source segment can be understood as a small segment of heat source in a homogeneous formation within the wellbore heat exchange section. Each heat source segment is located only within a single formation, therefore its thermal conductivity can be described using the equivalent thermophysical properties of the corresponding formation.

[0090] For example, a deep geothermal heat exchanger at a depth of 2500m passes through sandstone, limestone, and granite layers sequentially. After stratification in step S2, three finite linear heat source segments can be formed. The first heat source segment corresponds to the sandstone layer, the second to the limestone layer, and the third to the granite layer. Each heat source segment is described using the equivalent thermal conductivity and equivalent thermal diffusivity of its respective layer, without considering the influence of other strata. Subsequently, for each finite linear heat source segment, the corresponding transient thermal response is calculated under different time lag conditions. The transient thermal response refers to the temperature rise effect on any calculated node in the wellbore after a unit heat flow is released from a heat source segment and diffuses over a certain period of time. Since underground heat diffusion is a typical unsteady heat conduction process, the same heat source segment corresponds to different temperature rise responses at different times. For example, if heat is released only 1 hour later, the heat is mainly concentrated near the well wall; after 30 days, the heat has diffused to a more distant area, and its influence on the well wall temperature is significantly reduced. Therefore, it is necessary to calculate the thermal response for each time lag separately.

[0091] S312. Integrate and superimpose the transient thermal responses corresponding to each sub-path along the finite-line heat source integration path to obtain the thermal response matrix elements under the corresponding time lag between each calculation node and each heat source segment in the wellbore; wherein, for any calculation node With heat source section The thermal response matrix elements include real heat source items and surface-mirror heat source items, and the distance from the real heat source to the computing node. and the distance from the surface-mirror heat source to the computing node They are respectively: ; ; in, Where is the radius of the well wall. To calculate the node depth, For the heat source segment, the integral variable is... For any computing node With heat source section In time lag Under the condition, the corresponding thermal response matrix elements for: ; in, , These represent the actual heat source path at the th... The calculation point, the first Each heat source segment and Layered equivalent thermal conductivity and layered equivalent thermal diffusivity at each time step; , These represent the mirrored heat source paths at the [number]th ... The calculation point, the first Each heat source segment and Layered equivalent thermal conductivity and layered equivalent thermal diffusivity at each time step; Indicates a time lag. ; Indicates the time step; After completing the transient thermal response calculations for each heat source segment, step S312 further calculates the corresponding time-lag thermal response matrix elements between any calculation node in the wellbore and each heat source segment. Since the finite-line heat source model uses an integral form to describe the entire heat source, integration along the length of the heat source is required to accumulate the temperature responses generated at each location on the heat source segment in order to obtain the total temperature rise impact of that heat source segment on a specific calculation node. Therefore, this application integrates and superimposes the transient thermal responses corresponding to each sub-path along the finite-line heat source integration path to obtain the corresponding time-lag thermal response matrix elements between each calculation node in the wellbore and each heat source segment. Specifically, for any calculation node... and any heat source segment Its thermal response matrix elements include both real heat source terms and surface mirror heat source terms. The real heat source here refers to the finite linear heat source corresponding to the actual heat released from the wellbore; the mirror heat source is a mathematically equivalent heat source introduced to satisfy the surface adiabatic boundary conditions. Since the surface cannot absorb heat indefinitely, if only the real heat source is considered, the calculated temperature field cannot satisfy the actual boundary conditions. Therefore, this application adopts the mirror heat source method widely used in finite linear heat source theory, setting up a mirror heat source above the surface corresponding to the real heat source. Through the combined action of the real and mirror heat sources, the calculation results satisfy the surface boundary conditions, thereby improving the calculation accuracy of the analytical model.

[0092] S313. The thermal response matrix elements corresponding to each sub-path are summed to obtain the thermal response matrix elements between each calculation node and each heat source segment in the wellbore under the corresponding time lag. Step S313 further sums the thermal response matrix elements corresponding to each sub-path in the finite linear heat source integration path according to the integration path direction, thereby obtaining the thermal response matrix elements corresponding to the complete heat source segment. Since actual wellbore heat sources typically span multiple integration units, the temperature rise contribution of any heat source segment to the calculation node should be formed by the cumulative effect of all sub-paths, rather than using only the calculation result of a single integration point. Therefore, by summing the sub-path responses item by item, a more accurate overall thermal response of the heat source segment to the corresponding calculation node can be obtained.

[0093] S314. Store the thermal response matrix elements corresponding to each calculation node and each heat source segment in the wellbore at different time lags to form the thermal response matrix. Step S314 stores the thermal response matrix elements corresponding to all time lags, all calculation nodes, and all heat source segments in a unified order to form the thermal response matrix. This thermal response matrix can be understood as a "response database" of the underground surrounding rock heat transfer process. Each element in the matrix represents the contribution of a unit temperature rise to a calculation node after a unit heat flow is released from a heat source segment under a specific time lag condition. Subsequent step S32 does not require solving the underground heat conduction equation again. It only needs to call the thermal response matrix of the corresponding time lag and perform convolution summation with the heat flow increments released at each historical time step to quickly obtain the current wellbore temperature. Therefore, this application transforms the complex underground transient heat transfer calculation into matrix query and convolution operation by pre-computing the thermal response matrix. This not only ensures the physical accuracy of the layered stratum heat transfer calculation, but also significantly improves the computational efficiency of long-term operation prediction and large-scale simulation data generation. It provides an efficient and reliable physical prior foundation for the training and operation status prediction of subsequent physical-aided machine learning models.

[0094] In some embodiments of this application, S32 specifically includes: Let the increment of well depth heat flow at the nth time step be denoted as . Its heat flux increment is expressed as: ; where the initial time satisfies ; Based on the initial geothermal field and the thermal response matrix, the wellbore temperature is updated by performing historical heat flow convolution recursion using the historical heat flow increments: ; This is the wellbore temperature vector for all well depth calculation nodes at the nth time step; For the initial geothermal field The initial wellbore temperature vector is calculated based on the well depth nodes, and each component corresponds to the initial wellbore temperature at each calculation node. For time lag The corresponding thermal response matrix; For the first The historical well depth heat flow increment corresponding to each time step; S33 specifically includes: The wellbore temperature is used as the boundary condition for the coupled heat transfer equation of the wellbore fluid. Combined with the thermal resistance network between the wellbore fluid, the pipe wall and the surrounding rock, the discrete heat transfer equations for the annular fluid and the central pipe fluid are established by implicit upwind discretization. Wherein, the annular fluid temperature at the nth time step and the jth depth calculation node is denoted as The fluid temperature in the central tube is denoted as The annular fluid and the central tube fluid respectively satisfy the corresponding discrete heat transfer equations; The discrete equations for the annular fluid are: ; The discretization equation for the fluid in the central tube is: ; in, and These represent the equivalent linear heat capacities of the annular flow channel and the central pipe, respectively. This represents the measured mass flow rate of the circulating fluid; Indicates the specific heat capacity of the circulating fluid; Indicates the depth step size; This represents the thermal resistance between the annular fluid and the fluid in the central tube. Indicates the thermal resistance of the wellbore; The wellbore temperature at the j-th computation node at the n-th time step is represented by the wellbore temperature vector. The j-th component.

[0095] Boundary conditions are satisfied: ; This indicates the temperature of the annular fluid at the wellhead; Indicates the first Inlet temperature at each time step; The bottom of the well satisfies: ; This indicates the fluid temperature at the bottom of the well; This represents the annular fluid temperature at the bottom of the well; where N represents the total number of calculation nodes in the well depth direction. The outlet temperature meets the following requirements: ; Indicates the first Exit temperature at each time step; This indicates the fluid temperature at the center pipe of the wellhead; The discrete heat transfer equations for the annular fluid and the discrete heat transfer equations for the central tube fluid are combined into a linear equation set: ;in, This is the wellbore fluid coupling coefficient matrix; ; To form the load vector after moving the known terms in the discrete heat transfer equations of the annular fluid and the central tube fluid to the right side of the equations, the known terms include at least the wellbore temperature, the inlet boundary conditions, and the annular fluid temperature and the central tube fluid temperature at each depth calculation node in the previous time step. Solve the linear equations to obtain the temperature distribution of the fluid along the wellbore.

[0096] Specifically, in some embodiments of this application, step S33 is used to further calculate the temperature distribution of the circulating fluid inside the wellbore at various depths and time steps after obtaining the wellbore temperature, thereby establishing the coupled heat transfer relationship between the wellbore fluid and the surrounding rock, providing a basis for the subsequent calculation of operating parameters such as outlet temperature and heat exchange. For example, a deep geothermal heat exchanger can adopt a coaxial sleeve structure. After the circulating fluid enters the wellbore from the wellhead, it first flows towards the bottom of the well along the annulus channel. After reaching the bottom of the well, it turns back into the central tube and flows towards the wellhead along the central tube. Therefore, there are two flow paths inside the wellbore: the annulus fluid and the central tube fluid. During the fluid circulation process, the heat from the surrounding rock is first transferred to the annulus fluid through the wellbore wall. Subsequently, the annulus fluid and the central tube fluid continue to exchange heat through the central tube wall. At the same time, the fluid also continuously undergoes convective heat transfer along the axial direction. Therefore, the entire heat exchange process is a non-steady-state heat transfer process involving the coupling of the surrounding rock, the wellbore wall, the annulus fluid, and the central tube fluid.

[0097] For example, for a coaxial casing deep geothermal heat exchanger with a well depth of 3000m, the well depth can be discretized into multiple depth calculation nodes according to a fixed depth step, for example, discretized into 300 calculation nodes, each node corresponding to a well depth of 10m. For the j-th depth calculation node, at the n-th time step, the corresponding annular fluid temperature is established. and the temperature of the fluid in the central tube ,in, This represents the annular fluid temperature at the j-th depth calculation node at the n-th time step. This represents the center tube fluid temperature corresponding to the j-th depth calculation node at the nth time step. Subsequently, the wellbore temperature obtained in step S32 is used. As boundary conditions for the coupled heat transfer equations of the wellbore fluid, and in conjunction with the thermal resistance network formed between the wellbore, the central tube wall, and the fluid, discrete heat transfer equations corresponding to the annular fluid and the central tube fluid are established respectively.

[0098] The thermal resistance network is used to characterize the heat transfer path between the surrounding rock and the fluid inside the wellbore, and it mainly includes the wellbore thermal resistance. And the thermal resistance between the annular fluid and the central tube fluid. The wellbore thermal resistance This describes the thermal resistance generated during the transfer of heat from the surrounding rock through the wellbore to the annular fluid. Its magnitude can be calculated based on the thermal conductivity of the wellbore material, the wellbore thickness, and the heat transfer area. The thermal resistance is... It is used to describe the overall heat transfer resistance between the annular fluid and the central tube fluid. It typically includes the thermal resistance of the central tube wall and the convective heat transfer resistance of the fluid inside and outside the central tube. Therefore, it can comprehensively reflect the ease or difficulty of heat transfer from the annular fluid to the central tube fluid.

[0099] When establishing the fluid heat transfer equation, this application employs an implicit upwind discretization method to discretize and solve the fluid energy conservation equation. Implicit discretization means that all unknown temperatures within the current time step are treated as variables to be solved simultaneously in the equation, thus exhibiting good numerical stability and maintaining stable calculations even with large time steps. Upwind discretization means that the axial convection term is discretized using the upstream node temperature along the fluid flow direction, thereby conforming to the actual fluid flow characteristics. For example, when the annular fluid flows from the wellhead to the bottom, the annular fluid temperature corresponding to the j-th depth calculation node is mainly affected by the fluid temperature at the (j-1)-th depth calculation node; when the fluid in the central tube returns from the bottom to the wellhead, the discretization scheme is established using the upstream node temperature along the central tube flow direction, ensuring that the axial convection calculation conforms to the actual flow direction.

[0100] For annular fluids, the discrete heat transfer equation comprehensively considers four energy changes: the annular fluid's own heat capacity, heat transfer from the surrounding rock to the annular fluid, heat transfer from the annular fluid to the fluid in the central tube, and axial convection of the fluid. Among these, It represents the equivalent linear heat capacity corresponding to the annular flow channel, used to characterize the ability of the annular fluid to store heat per unit well depth; This indicates the mass flow rate of the circulating fluid, which can be obtained by real-time measurement using a flow meter. This indicates the specific heat capacity of the circulating fluid, which can be determined according to the type of circulating medium. For example, the specific heat capacity of water is approximately 4200 J / (kg·K). Indicates the depth step size in the well depth direction; This represents the wellbore temperature at the corresponding depth calculation node, which is calculated using the thermal response matrix in step S32. Therefore, the discrete heat transfer equation for the annular fluid can simultaneously reflect the continuous heating of the surrounding rock to the fluid, the fluid's own thermal inertia, and the heat migration process brought about by axial flow.

[0101] Similarly, for fluid in a central tube, its discrete heat transfer equation comprehensively considers factors such as the heat capacity of the fluid in the central tube itself, heat transfer from the annular fluid to the fluid in the central tube, and axial convection of the fluid inside the central tube. Among these factors, This represents the equivalent linear heat capacity corresponding to the central tube, used to characterize the heat storage capacity of the fluid per unit well depth inside the central tube. Since the fluid inside the central tube mainly receives heat transferred from the annular fluid, the temperature of the central tube fluid is affected not only by the temperature of the central tube fluid at the previous node, but also by the temperature of the annular fluid at the corresponding depth, thus forming a two-way coupling relationship between the annular fluid and the central tube fluid.

[0102] To ensure that the entire wellbore fluid heat transfer process meets actual operating conditions, this application further sets corresponding boundary conditions. Specifically, the wellhead inlet boundary satisfies... That is, the annular fluid temperature at the wellhead is equal to the inlet temperature corresponding to the current time step. For example, if the temperature of the circulating fluid entering the wellhead is 18°C, then... =18℃. The bottom connection boundary satisfies... This means that because the annular fluid deflects back into the central tube at the bottom of the well, the temperature of the fluid in the central tube at the bottom of the well remains continuous with the temperature of the fluid in the annular fluid at the bottom of the well. The outlet boundary satisfies... The fluid temperature in the central pipe at the wellhead is the final outlet temperature of the deep geothermal heat exchanger, which is also an important output parameter for predicting subsequent operating conditions.

[0103] After establishing the discrete equations for the annular fluid and the central tube fluid corresponding to all depth calculation nodes, all discrete equations can be unified into a linear equation system: .

[0104] in, This represents the wellbore fluid coupling coefficient matrix, which consists of the heat transfer coefficient and wellbore thermal resistance corresponding to each calculation node. Thermal resistance of heat transfer between fluids Circulating fluid mass flow rate Specific heat capacity of circulating fluid Equivalent linear heat capacity and Depth step size and time step The matrix is ​​determined jointly, thus it can comprehensively reflect the coupling relationship between nodes at various depths inside the wellbore and between the annular fluid and the central tube fluid.

[0105] The unknown vector It consists of all the fluid temperatures to be determined at the current time step, that is: .

[0106] For example, when the well depth direction is discretized into 300 computational nodes, the unknown vector It contains a total of 602 unknown temperature variables, including 301 annular fluid temperatures and 301 central tube fluid temperatures.

[0107] The load vector To unify and move all known terms in the discrete heat transfer equations of the annular fluid and the discrete heat transfer equation of the central tube fluid to the right-hand side of the equation, the resulting vector mainly consists of three parts: the first part is the wellbore temperature corresponding to each depth calculation node at the current time step, obtained from step S32. The resulting surrounding rock heating term represents the heat continuously input from the surrounding rock to the wellbore fluid; the second part is the wellhead inlet temperature. Boundary constraints, such as the continuity of bottom hole temperature, are used to ensure that the wellbore fluid flow meets the actual operating conditions; the third part is the annular fluid temperature corresponding to each depth calculation node in the previous time step. and the temperature of the fluid in the central tube The resulting historical temperature term is used to reflect the fluid's thermal inertia and heat accumulation effect. For example, if the annular fluid temperature at a certain depth in the previous hour was 35°C, this temperature is used as a known term in the calculation of the corresponding node temperature for the next time step, making the fluid temperature change more consistent with the actual continuous heat transfer law. Therefore, the load vector... This comprehensively reflects the combined effects of the current surrounding rock heating effect, inlet operating conditions, and historical thermal state on the temperature distribution of the wellbore fluid.

[0108] Finally, Gaussian elimination, LU decomposition, chasing method, or sparse matrix direct solution algorithm can be used to solve the linear equation system. Numerical solutions are performed to obtain the annular fluid temperature and central tube fluid temperature at each calculation node along the depth direction of the wellbore at the current time step, thus obtaining the complete temperature distribution of the wellbore fluid along the friction path. Since the entire solution process simultaneously considers multiple influencing factors such as surrounding rock heat transfer, wellbore wall heat conduction, heat exchange between fluids, axial convection of the fluid, and historical thermal inertia, the obtained temperature distribution of the wellbore fluid along the friction path can more realistically reflect the actual heat transfer state of the deep geothermal heat exchanger under long-term operating conditions. This provides a reliable basis for calculating the output parameters of the analytical heat transfer model, such as outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term thermal decay index.

[0109] Figure 2This is a schematic diagram of the wellbore structure and layered formations applicable to deep coaxial casing geothermal heat exchangers according to the present invention. Figure 2 As shown, the heat exchange well traverses multiple strata with different thermal properties along its depth, including soil geological layer 1, soil geological layer 2, soil geological layer 3, and soil geological layer 4. Each layer has different thermal conductivity, density, specific heat capacity, and geothermal gradient. The well depth is typically greater than 2000 m, and the initial underground temperature gradually increases along the well depth according to the geothermal gradient. Figure 2 middle This represents the geothermal gradient. The wellbore employs a coaxial casing structure, consisting of a central tube, annular channel, outer tube, and wellbore backfill. During operation, the cryogenic working fluid enters the annular channel through the wellhead inlet, flows downwards along the wellbore, reverses direction at the bottom of the well, enters the central tube, then flows upwards along the central tube, and exits at the wellhead at the outlet temperature. During this flow, the working fluid exchanges heat with the surrounding layered rock through the outer tube and backfill material. The thermophysical parameters of strata at different depths collectively affect the heat transfer characteristics and heat exchange efficiency of the wellbore, thus forming a complex wellbore-surrounding rock coupled heat transfer process. This diagram illustrates... Figure 2 This reflects the layered formation, coaxial wellbore, and bidirectional flow heat transfer model established in the embodiments of this application, which provides a physical model basis for subsequent layered finite line heat source response calculation, thermal resistance network establishment, and wellbore fluid temperature solution.

[0110] Example 2; The analytical heat transfer model established in this application is a layered coupled rapid analytical model applicable to deep coaxial casing geothermal heat exchangers. The analytical heat transfer model uses wellbore structural parameters, drilling exploration data, and operating condition parameters as inputs. The model inputs include at least: well depth. Inner and outer pipe diameters , , , (in, Indicates the inner diameter of the inner tube (i.e., the central tube). Indicates the outer diameter of the inner tube; Indicates the inner diameter of the outer annulus pipe. Indicates the outer diameter of the outer tube and the diameter of the borehole. Thermal conductivity of pipe wall and backfill material, fluid density Specific heat capacity Dynamic viscosity inlet temperature at time t Volumetric flow rate Operating heat load Time step Surface temperature and the top and bottom depths, thermal conductivity, density, specific heat capacity, and first degree of heat of each stratum. Geothermal gradient of the layer The model output includes the outlet temperature at time t. Heat transfer per unit well depth, cumulative heat transfer, wellbore temperature, and long-term thermal decay index .

[0111] With depth Downward is the positive direction, if the first The top and bottom depths of the strata are and And geothermal gradient Expressed in °C / m, the original geothermal field is defined as follows:

[0112] in, For depth The initial ground temperature at that location, This is the surface reference temperature. For the first Layer thickness, For the first Geothermal gradient, For the first Geothermal gradient; Temperature continuity is satisfied at the interfaces of each layer: , Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time; Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time. The heat flow continuity condition is satisfied at each level: ,in, Indicates the first The thermal conductivity of the strata, Indicates the first The thermal conductivity of the strata, Indicates from the first One side of the layer approaches the interface Temperature gradient over time; Indicates from the first One side of the layer approaches the interface The temperature gradient over time; n is the normal direction of the interlayer interface; Indicates from the first One side of the layer approaches the interlayer interface. Indicates from the first One side of the layer approaches the interlayer interface.

[0113] For heat transfer paths crossing one or more formation interfaces, this application embodiment does not directly use a single homogeneous formation parameter, but instead divides the path into several sub-paths according to the formation interfaces, that is, for the equivalent thermal conductivity... equivalent density Equivalent specific heat capacity and equivalent thermal diffusivity The calculation method corresponds to the simplification of interlayer coupling boundaries when a finite line heat source traverses a layered medium, and can preserve the influence of differences in the thermal properties of different strata on the heat transfer path. This improves the model's applicability to complex layered strata and its computational accuracy.

[0114] A thermal resistance network is used to couple the wellbore fluid to the pipe wall. The equivalent thermal resistance between the annular fluid and the central pipe fluid is... And the equivalent thermal resistance between the annular fluid and the surrounding rock outside the wellbore. , respectively represented as: ; ; in, and These are the convective heat transfer coefficients on the central tube side and the annulus side, respectively. , and These are the thermal conductivity coefficients of the central pipe wall, outer pipe wall, and backfill layer, respectively. The convective heat transfer coefficient can be updated in real time by the flow rate. Volumetric flow rate is converted to mass flow rate. The flow velocity, Reynolds number Re, and Nusselt number Nu are calculated based on the flow area of ​​the central tube and the flow area of ​​the annulus. Under turbulent conditions, the Gnielinski relation is used. Taking under laminar flow conditions , This indicates the friction factor of the flow inside the pipe. .Depend on Get different traffic and This allows the thermal resistance network to be dynamically updated according to operating conditions.

[0115] Unsteady-state heat transfer in soil is studied using a stratified finite-line heat source response. The thermal response of the surrounding rock is divided along the well depth direction into... Each heat source segment, calculation point Located at depth heat source section The inner integral variable is The radius of the well wall is Then time lag The response matrix elements below It can be defined as: ; The first item is the real heat source item, and the second item is the surface-image heat source item. , These represent the actual heat source path at the th... The calculation point, the first The heat source segment and the first Layered equivalent thermal conductivity and layered equivalent thermal diffusivity under time lag; , These represent the mirrored heat source paths at the [number]th ... The calculation point, the first The heat source segment and the first Layered equivalent thermal conductivity and layered equivalent thermal diffusivity under time lag; Indicates a time lag. ; Indicates the time step.

[0116] During implementation, for each heat source segment and calculation points First, calculate and store the sub-path length using Gaussian integral points. , and their corresponding equivalent thermophysical parameters, then according to time lag Pre-form a set of response matrices Therefore, during long-term simulations, the ERFC integral is no longer repeatedly calculated; instead, the pre-calculated matrix is ​​used to complete the thermal response superposition.

[0117] The historical heat flow convolution recursion is specifically as follows: Let the first... The heat flux vector per unit well depth at each time step is , No. The heat flux increment at each time step is ,and Then the wellbore temperature vector is updated to Each time step forward will bring the latest... The heat flux history queue is pushed in, and the response matrix is ​​multiplied and summed with the heat flux increment queue in reverse time lag order to obtain the current wellbore temperature. This recursive method directly corresponds to the principle of superposition of heat conduction and avoids re-integrating all historical heat fluxes.

[0118] In obtaining wellbore temperature Then, the implicit upwind scheme is used to solve for the wellbore fluid temperature. Let... and The first The time step, the first Annular fluid temperature and central tube fluid temperature at each depth node and These are the equivalent linear heat capacities of the annular channel and the central tube channel, respectively. Given a depth step size, the discretization equations for the annular fluid and the central tube fluid are respectively: ; Boundary conditions are Bottom of the well The outlet temperature is The above system of equations can be written in matrix form. ,in The temperature distribution along the path at the current moment is obtained by solving the matrix. Among them, This indicates the temperature of the annular fluid at the wellhead; Indicates the first The outlet temperature at each time step; This indicates the fluid temperature at the bottom of the well; This represents the annular fluid temperature at the bottom of the well; where N represents the total number of calculation nodes in the well depth direction. Indicates the first The outlet temperature at each time step; This represents the fluid temperature at the center pipe of the wellhead. The cumulative heat transfer is obtained by time integration. After the above calculations are completed, update... Then proceed to the next time step, forming a closed-loop rapid calculation process of "response matrix pre-calculation - wellbore temperature convolution recursion - wellbore fluid matrix solution - heat flux increment update".

[0119] For parameterized sample generation and physics-assisted machine learning training, since the number of actual test samples is usually small and the coverage of working conditions is limited, this application will use simulated sample sets. Compared with the measured sample set These are used together for training. First, the range of values ​​for each variable is determined, and initial candidate working conditions are obtained using Latin hypercube, Sobol low-difference sequences, and equally spaced stratified sampling. Then, a score is constructed based on the rate of change output by the physical model, the distance between candidate samples, the prediction uncertainty of the surrogate model, and whether it belongs to an extreme working condition. .

[0120] The simulated dataset generated by the above intelligent sampling strategy can be represented as follows: ,in This is the system and runtime input vector (where the subscript i represents the sample number). For the parsing model output, To exchange heat, As a long-term thermal decay index, Used to label regular, boundary, and extreme sample types.

[0121] In one specific implementation, the data-driven model employs a physical residual encoder-multiple decoder architecture, rather than a single fully connected prediction network. Static encoder Used to encode well depth, pipe diameter, thermal resistance, formation thickness, thermal conductivity, density, specific heat capacity, geothermal gradient, and statistical characteristics of the response matrix; time-series encoder Used to encode time series such as inlet temperature, flow rate, heat load, start / stop status, and historical outlet temperature; physical prior branch. The input consists of the predicted sequence and response kernel features of the analytical model. The above encoding results are then fused to form latent variables. ,Right now This is represented as the fused latent feature vector.

[0122] Multiple decoders include an outlet temperature decoder. Heat exchanger decoder Long-term thermal attenuation decoder And measured deviation correction decoder The final output of the model is represented as The prediction of medium-deep geothermal heat exchangers includes both static formation heterogeneity and strong temporal operating boundaries. The encoder can compress high-dimensional layered parameters, thermal resistance parameters, and operating history into a unified state, while the decoder can output time series or engineering indicators for different physical quantities, thus matching the application scenarios of outlet temperature, heat exchange, and long-term decay.

[0123] The training objective function can be expressed as: .in To simulate sample residual error, This is the error of the measured sample. To smooth out errors over time, Error correction is applied to the deviation between the simulated and measured domains. Time smoothing constraints are used to suppress physical oscillations under non-extreme conditions; measured sample weights... It can be higher than the simulated sample weights. Alternatively, adjustments can be made dynamically based on sample noise estimation.

[0124] Regarding model evaluation and subsequent online updates: Using measured outlet temperature, measured heat transfer, and distributed temperature test results, the physical enhancement data-driven model output is compared with field data, employing root mean square error (RMSE). and the coefficient of determination The evaluation is conducted, and the calculation methods are as follows: ; ; Where N is the total amount of data in the dataset. For predicted values, These are measured values. This represents the average value in the dataset.

[0125] When a new operating condition arrives, a joint prediction is first performed using a fast analytical model and a physical residual learning model. If the evaluation index under the new operating condition exceeds a threshold, local perturbation samples are generated centered on that operating condition, and the bias correction decoder or end fusion layer is updated first. When the cumulative number of new samples reaches a set number, the static encoder and the temporal encoder are then jointly fine-tuned with a small learning rate. This enables the model to continuously absorb long-term field operation data while reducing the risk of overfitting caused by frequent retraining of the entire model.

[0126] The method described in this application can rapidly obtain the outlet temperature, instantaneous heat transfer, cumulative heat recovery, wellbore temperature changes, and long-term thermal decay trends under complex stratified formations, small sample measurements, and long-term operation conditions. This method is applicable not only to heat exchanger selection and design, parameter optimization, and operation control, but also to engineering project scheme comparison, life assessment, digital twin platform construction, and operational anomaly early warning. This invention is not limited to coaxial sleeve structures; it can also be applied to U-tube or other deep geothermal heat exchanger structures by adjusting the fluid equations, thermal resistance networks, and response matrix construction.

[0127] Although embodiments of this application have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting this application. Those skilled in the art can make modifications, alterations, substitutions and variations to the above embodiments within the scope of this application.

Claims

1. A method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning, characterized in that, include: S1. Obtain the wellbore structural parameters, drilling exploration data, and operating condition parameters of the deep geothermal heat exchanger; establish an analytical heat transfer model considering soil and rock stratification and geothermal gradient based on the drilling exploration data, and determine the initial geothermal field of the analytical heat transfer model based on the geothermal gradient of each stratum. S2. Based on the analytical heat transfer model, the finite linear heat source integral path is processed in layers, and the corresponding equivalent thermal properties are calculated according to the thermal properties of each layer, thereby establishing a layered thermal response relationship that satisfies the interlayer temperature continuity condition and the heat flow continuity condition. S3. Calculate the thermal response matrix based on the layered thermal response relationship, and solve the wellbore temperature by combining the initial geothermal field and historical heat flow increment; establish the coupled heat transfer relationship between the wellbore fluid and the surrounding rock based on the wellbore temperature, and obtain the output of the analytical heat transfer model; S4. Using the analytical heat transfer model, parameter calculations are performed under different wellbore structure parameters, operating condition parameters, and drilling exploration data to generate simulated operating data. Based on the simulated operating data, adaptive enhancement is performed to construct a simulated dataset. S5. Use the simulated dataset and the measured running data to train the physics-assisted machine learning model and obtain the trained physics-assisted machine learning model. S6. Obtain the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted. Based on the wellbore structural parameters, drilling exploration data, and operating condition parameters corresponding to the operating condition to be predicted, obtain the output of the analytical heat transfer model corresponding to the operating condition to be predicted. Input the wellbore structural parameters, drilling exploration data, operating condition parameters, and the output of the analytical heat transfer model into the trained physical-aided machine learning model to obtain the predicted operating status of the deep geothermal heat exchanger.

2. The method for predicting the operating status of deep geothermal heat exchangers using physical-aided machine learning according to claim 1, characterized in that, The method further includes: S7. Acquire new measured operating data generated during the operation of the deep geothermal heat exchanger, and use the trained physical-aided machine learning model to predict the operating conditions corresponding to the new measured operating data. When the prediction error between the prediction result and the measured result corresponding to the newly added measured operating data exceeds a preset threshold, new simulated operating data is generated using the operating conditions corresponding to the newly added measured operating data, and the new simulated operating data and the new measured operating data are used to incrementally update the physics-assisted machine learning model.

3. The method for predicting the operating status of deep geothermal heat exchangers using physical-aided machine learning according to claim 1, characterized in that, The wellbore structural parameters include at least the well depth, well diameter, well wall radius, wellbore heat exchange section length, central tube inner diameter, central tube outer diameter, annulus inner diameter, annulus outer diameter, casing inner diameter, casing outer diameter, central tube wall thickness, and casing wall thickness. The drilling exploration data includes at least the top and bottom depths, layer thicknesses, thermal conductivity, density, specific heat capacity, thermal diffusivity, and geothermal gradients of each stratum; The operating parameters include at least inlet temperature, circulation flow rate, heat load, and operating time; The initial geothermal field satisfies: ； in, For depth The initial ground temperature at that location, This is the surface reference temperature. For the first Geothermal gradient of the strata For the first Geothermal gradient of the strata For the first The thickness of the strata, and The first The depth of the top and bottom interfaces of the strata; The initial geothermal field satisfies the temperature continuity condition at each stratum interface: ; Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time; Indicates the first One side of the layer approaches the interlayer interface The initial ground temperature at that time; where the superscript "-" indicates approaching the layer from the top side of the interface, and the superscript "+" indicates approaching the interface from the bottom side of the interface; And satisfy the heat flow continuity condition: ; in, Indicates the first The thermal conductivity of the strata Indicates the first The thermal conductivity of the strata Indicates from the first One side of the layer approaches the interface Temperature gradient over time; Indicates from the first One side of the layer approaches the interface The temperature gradient over time; n is the normal direction of the interlayer interface; Indicates from the first One side of the layer approaches the interlayer interface. Indicates from the first One side of the layer approaches the interlayer interface.

4. The method for predicting the operating status of deep geothermal heat exchangers using physical-aided machine learning according to claim 3, characterized in that, S2 specifically includes: S21. Determine the finite linear heat source integration path based on the positional relationship between the wellbore heat exchange section and the layered formation; S22. When the finite linear heat source integration path crosses a formation interface, the finite linear heat source integration path is segmented according to the formation interface to form multiple sub-paths located in different formations. S23. Calculate the equivalent thermal properties of each sub-path based on the thermal properties of the formation corresponding to each sub-path; The thermal properties include thermal conductivity, density, and specific heat capacity; the equivalent thermal properties include: equivalent thermal conductivity. equivalent density Equivalent specific heat capacity and equivalent thermal diffusivity ;in, ; ; ; ; To indicate the first The proportion of each sub-path length to the integral path length of the finite line heat source. , and They represent the first The thermal conductivity, density, and specific heat capacity of the strata corresponding to each sub-path; Indicates the subpath number involved in the summation; Indicates the first The proportion of the length of each sub-path to the length of the integral path of the finite linear heat source; Indicates the first The density of the stratum where each sub-path is located; S24. Calculate the temperature response function of the finite linear heat source in the corresponding sub-path based on the equivalent thermal property parameters of each sub-path. Combine the temperature response functions according to the position order of each sub-path in the integral path of the finite linear heat source, and satisfy the temperature continuity condition and heat flow continuity condition of the adjacent strata interface to establish the layered thermal response relationship of the finite linear heat source in the layered strata.

5. The method for predicting the operating status of deep geothermal heat exchangers using physical-aided machine learning according to claim 4, characterized in that, Step S3 specifically includes: S31. Based on the hierarchical thermal response relationship, the transient heat transfer response of the finite linear heat source under different time lags is pre-calculated, the thermal response function between each calculation node of the wellbore and each heat source segment is calculated, and the thermal response matrix is ​​constructed according to the thermal response function corresponding to each time lag. S32. Based on the initial geothermal field, the thermal response matrix, and the historical heat flux increment, calculate the wellbore temperature corresponding to each time step by recursive calculation through historical heat flux convolution. S33. Using the wellbore temperature as the boundary condition for the coupled heat transfer equation of the wellbore fluid, and combining it with the thermal resistance network between the wellbore fluid, the pipe wall and the surrounding rock, solve for the temperature distribution of the wellbore fluid along the pipe. S34. Determine the output of the analytical heat transfer model based on the temperature distribution along the wellbore fluid. The output of the analytical heat transfer model includes at least the outlet temperature, heat transfer per unit well depth, instantaneous heat transfer, and long-term heat decay index.

6. The method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning according to claim 5, characterized in that, S31 specifically includes: S311. Based on the layered thermal response relationship, the multiple sub-paths corresponding to the finite linear heat source integral path are respectively taken as finite linear heat source segments in the corresponding strata, and the transient thermal response under the corresponding time lag is calculated using the equivalent thermal conductivity and equivalent thermal diffusivity corresponding to each sub-path. S312. Integrate and superimpose the transient thermal responses corresponding to each sub-path along the finite-line heat source integration path to obtain the thermal response matrix elements under the corresponding time lag between each calculation node and each heat source segment in the wellbore; wherein, for any calculation node With heat source section The thermal response matrix elements include real heat source items and surface-mirror heat source items, and the distance from the real heat source to the computing node. and the distance from the surface-mirror heat source to the computing node They are respectively: ； ； in, Where is the radius of the well wall. To calculate the node depth, For the integral variable of the heat source segment; For any computing node With heat source section In time lag Under the condition, the corresponding thermal response matrix elements for: ； in, , These represent the actual heat source path at the th... The calculation point, the first Each heat source segment and Layered equivalent thermal conductivity and layered equivalent thermal diffusivity at each time step; , These represent the mirrored heat source paths at the [number]th ... The calculation point, the first Each heat source segment and Layered equivalent thermal conductivity and layered equivalent thermal diffusivity at each time step; Indicates a time lag. ; Indicates the time step; S313. Accumulate the thermal response matrix elements corresponding to each sub-path to obtain the thermal response matrix elements between each calculation node and each heat source segment in the wellbore with the corresponding time lag. S314. Store the thermal response matrix elements corresponding to each calculation node and each heat source segment in the wellbore with different time lags to form the thermal response matrix.

7. The method for predicting the operating status of deep geothermal heat exchangers using physical-aided machine learning according to claim 6, characterized in that, Specifically, S32 includes: Let the increment of well depth heat flow at the nth time step be denoted as... Its heat flux increment is expressed as: ; where the initial time satisfies ; Based on the initial geothermal field and the thermal response matrix, the wellbore temperature is updated by performing historical heat flow convolution recursion using the historical heat flow increments: ; This is the wellbore temperature vector for all well depth calculation nodes at the nth time step; For the initial geothermal field The initial wellbore temperature vector is calculated based on the well depth nodes, and each component corresponds to the initial wellbore temperature at each calculation node. For time lag The corresponding thermal response matrix; For the first The historical well depth heat flow increment corresponding to each time step; S33 specifically includes: The wellbore temperature is used as the boundary condition for the coupled heat transfer equation of the wellbore fluid. Combined with the thermal resistance network between the wellbore fluid, the pipe wall and the surrounding rock, the discrete heat transfer equations for the annular fluid and the central pipe fluid are established by implicit upwind discretization. Wherein, the annular fluid temperature at the nth time step and the jth depth calculation node is denoted as The fluid temperature in the central tube is denoted as The annular fluid and the central tube fluid respectively satisfy the corresponding discrete heat transfer equations; The discrete equations for the annular fluid are: ； The discretization equation for the fluid in the central tube is: ； in, and These represent the equivalent linear heat capacities of the annular flow channel and the central pipe, respectively. This represents the measured mass flow rate of the circulating fluid; Indicates the specific heat capacity of the circulating fluid; Indicates the depth step size; This represents the thermal resistance between the annular fluid and the fluid in the central tube. Indicates the thermal resistance of the wellbore; The wellbore temperature at the j-th computation node at the n-th time step is represented by the wellbore temperature vector. The j-th component; Boundary conditions are satisfied: ; This indicates the temperature of the annular fluid at the wellhead; Indicates the first Inlet temperature at each time step; The bottom of the well satisfies: ; This indicates the fluid temperature at the bottom of the well; This represents the annular fluid temperature at the bottom of the well; where N represents the total number of calculation nodes in the well depth direction. The outlet temperature meets the following requirements: ; Indicates the first Exit temperature at each time step; This indicates the fluid temperature at the center pipe of the wellhead; The discrete heat transfer equations for the annular fluid and the discrete heat transfer equations for the central tube fluid are combined into a linear equation set: ;in, This is the wellbore fluid coupling coefficient matrix; ； To form the load vector after moving the known terms in the discrete heat transfer equations of the annular fluid and the discrete heat transfer equation of the central tube fluid to the right side of the equations, the known terms include at least the wellbore temperature, the inlet boundary conditions, and the annular fluid temperature and the central tube fluid temperature at each depth calculation node in the previous time step. This represents a column vector consisting of all the temperatures to be solved at the current time step; Solve the linear equations to obtain the temperature distribution of the fluid along the wellbore.

8. The method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning according to claim 7, characterized in that, S4 specifically includes: S41. Based on the wellbore structure parameters, operating condition parameters and drilling exploration data, determine the input parameters and the value range of each input parameter for the analytical heat transfer model. The initial simulation conditions are generated by combining Latin hypercube sampling, Sobol low-difference sequence sampling and equally spaced layered sampling. S42. Input the initial simulation conditions into the analytical heat transfer model to obtain the analytical heat transfer model output corresponding to each initial simulation condition, and generate the corresponding initial simulation operation data. S43. Based on the output of the analytical heat transfer model, determine the physical response sensitivity of each input parameter of the analytical heat transfer model, and construct a sampling scoring function by combining the prediction uncertainty of the surrogate model and the boundary distance of the operating condition. Based on the sampling scoring function, perform adaptive sampling on the parameter space to obtain the new simulated operating condition. The sampling scoring function is expressed as follows: ； in, This indicates the uncertainty in the surrogate model prediction. This represents the physical response sensitivity corresponding to the output of the analytical heat transfer model. Indicates the distance to the working condition boundary. Indicates indicators under extreme operating conditions. These represent the corresponding weights; For sampling scoring functions; S44. Input the newly added simulated operating conditions and the preset extreme operating conditions into the analytical heat transfer model to generate corresponding newly added simulated operating data, and combine the newly added simulated operating data with the initial simulated operating data to form the simulation dataset.

9. The method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning according to claim 8, characterized in that, The physics-assisted machine learning model includes a static encoder, a temporal encoder, a physical prior encoder, and multiple decoders, wherein: The static encoder is used to encode wellbore structural parameters and drilling exploration data; The timing encoder is used to encode timing data corresponding to inlet temperature, flow rate, heat load, start / stop status, and historical operating status. The physical prior encoder is used to encode the output of the parsed heat transfer model and the response kernel features extracted from the thermal response matrix; The features output by the static encoder, the temporal encoder, and the physical prior encoder are fused to form a unified latent variable: ;in, Represents a unified latent variable. This represents the static coding features output by the static encoder. This represents the timing coding characteristics output by the timing encoder. F represents the physical prior encoded features output by the physical prior encoder, and F() represents the feature fusion function; The unified latent variable Multiple decoders are input, and the outlet temperature, heat exchange, long-term heat decay index and prediction deviation correction amount are output respectively to obtain the corrected deep geothermal heat exchanger operation prediction model. The multiple decoders include an outlet temperature decoder, a heat exchanger decoder, a long-term thermal decay decoder, and a deviation correction decoder. The outlet temperature decoder is used to output the outlet temperature of the deep geothermal heat exchanger. The heat exchanger decoder is used to output the heat exchanger of the deep geothermal heat exchanger. The long-term thermal decay decoder is used to output the long-term thermal decay index of the deep geothermal heat exchanger. The deviation correction decoder is used to learn and analyze the prediction deviation between the output of the heat transfer model and the measured operation results, and output the corresponding prediction deviation correction amount. The predicted deviation correction amount is fused with the output of the analytical heat transfer model to obtain the corrected prediction result for the deep geothermal heat exchanger operation, which is expressed as follows: ; in, This represents the output of the analytical heat transfer model. This indicates a deviation correction decoder. The corrected operating prediction results for the deep geothermal heat exchanger include outlet temperature, heat exchange capacity, and long-term heat decay index.

10. The method for predicting the operating status of a deep geothermal heat exchanger using physical-aided machine learning according to claim 9, characterized in that, The training objective function of the physics-assisted machine learning model is expressed as: ; in, This represents the prediction error corresponding to the simulated dataset; This represents the prediction error corresponding to the measured operating data; This represents the time series smoothing constraint error of the time series data corresponding to the operating condition parameters; This represents the domain deviation correction error between the simulated dataset and the measured operating data; , , and These represent the weighting coefficients of the corresponding loss terms; For training objective function; When the training objective function converges to a preset loss threshold, or when the change in the training objective function within a preset number of consecutive training rounds is less than a preset convergence threshold, the training of the physical-assisted machine learning model is determined to be complete, and a trained physical-assisted machine learning model is obtained.

Images