A Dynamic Heating Method and System for Frost Heave Prevention in Cold Region Tunnels Based on Multiphysics Field Coupling

Abstract

This invention provides a dynamic heating method and system for preventing frost heave in cold-region tunnels based on multi-physics coupling, relating to the field of frost heave prevention technology for cold-region tunnels. The method includes acquiring temperature gradients, moisture content distributions, and freezing front advance speeds at different depths within the tunnel; performing coupling analysis to identify critical frost heave risk zones; determining heating response sequences; identifying heat transfer path networks; calculating heat utilization efficiency coefficients; and generating dynamic heating source spatial configuration schemes. This invention, through multi-physics coupling analysis and a dynamic response mechanism, accurately locates heating positions and timing, improves energy utilization efficiency, and effectively prevents frost heave damage.

Description

Technical Field

[0001] This invention relates to the field of frost heave prevention technology for tunnels in cold regions, and particularly to a dynamic heating method and system for frost heave prevention in tunnels in cold regions based on multi-physics field coupling. Background Technology

[0002] In cold-region tunnel engineering, frost heave is a severe challenge. When temperatures drop below freezing, moisture in the soil or rock freezes and expands, generating frost heave force that severely damages the tunnel structure. Tunnels in cold regions affected by frost heave often suffer from lining cracking, deformation, and leakage, seriously threatening the structural safety and operational stability of the tunnel. Currently, frost heave prevention technologies for cold-region tunnels mainly include thermal insulation, electric heating, and geothermal water circulation. Among these, electric heating has become an important means of preventing frost heave in cold-region tunnels due to its advantages such as flexible control and convenient implementation. Traditional electric heating systems mainly use heating wires and cables, providing heat to resist frost damage by pre-embedding heating elements in the lining or surrounding rock.

[0003] Traditional heating systems generally employ constant temperature control or simple temperature threshold triggering modes, failing to fully consider the coupling effect of moisture migration and temperature field in the tunnel surrounding rock. This results in ineffective anti-freezing measures and a lack of specificity in addressing dynamic changes in areas at risk of frost heave. Existing heating solutions mostly rely on experience to determine heating locations and power configurations, lacking systematic analysis of heat transfer paths. This hinders efficient heat utilization, leading to energy waste and unsatisfactory heating results. Traditional anti-frost heave systems fail to establish a temporal correlation between the dynamic advancement of freezing fronts and heating response, making precise control of heating measures difficult. Under the complex coupling effects of temperature and moisture fields, they cannot respond promptly and effectively to changes in frost heave risk in different areas, compromising the effectiveness of anti-freezing measures. Summary of the Invention

[0004] This invention provides a dynamic heating method and system for preventing frost heave in cold-region tunnels based on multi-physics field coupling, which can solve the problems in the prior art.

[0005] A first aspect of the present invention provides a dynamic heating method for preventing frost heave in cold-region tunnels based on multi-physics field coupling, comprising:

[0006] To obtain the temperature gradient, water content distribution, and freezing front advance velocity at different depths in tunnels in cold regions;

[0007] The coupling analysis of the temperature gradient and the moisture content distribution is used to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and the distribution map of the critical area of ​​frost heave risk is obtained.

[0008] Based on the freezing front advance speed of each segment in the distribution map of critical areas of freezing heave risk, a time priority ranking is established by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface, and a heating response sequence with time dimension constraints is obtained.

[0009] Based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracking the evolution trajectory of the temperature isosurface, and the heat transfer path network topology is obtained.

[0010] In each path of the transmission path network topology, the heat utilization efficiency coefficient of different heating positions in each section is calculated by quantifying the contribution of temperature rise and energy attenuation when a unit heating power is transmitted to the freezing front along the current path.

[0011] Based on the heating response sequence and the heat utilization efficiency coefficient, a dynamic heating source spatial configuration scheme for frost heave prevention in cold region tunnels is generated.

[0012] The coupling analysis of the temperature gradient and the moisture content distribution yields the driving potential field for moisture migration. Based on the spatial distribution characteristics of this driving potential field, risk zones in cold-region tunnels where moisture continuously migrates and accumulates towards the freezing front are identified, resulting in a distribution map of critical areas for frost heave risk, including:

[0013] Establish a spatial discrete grid in the three-dimensional space of tunnels in cold regions;

[0014] The three-dimensional spatial components of the temperature gradient and the numerical distribution of the water content are registered and aligned on the same spatial discrete grid nodes. The water migration driving potential field is obtained by calculating the superposition effect of the adsorption potential of ice-water phase transition caused by temperature decrease and the osmotic potential generated by the water content gradient at each grid node.

[0015] The spatial gradient of potential energy value is extracted from the water migration driving potential field. The direction of the spatial gradient is taken as the preferred water migration direction, and the magnitude of the spatial gradient is taken as the water migration driving intensity, so as to obtain a potential field distribution vector field containing direction information and intensity information.

[0016] In the potential field distribution vector field, trace the potential energy gradient streamlines along the direction of potential energy reduction from the region with water content greater than the preset reference value to the region with temperature below the freezing point, and identify the convergence region where multiple potential energy gradient streamlines intersect in space.

[0017] Calculate the convergence density and convergence intensity of the potential energy gradient streamlines within the convergence region, and mark the convergence region where the convergence density exceeds a preset density threshold as a risk zone for continuous water migration and accumulation.

[0018] The location coordinates, spatial range, and convergence intensity of the risk section in the three-dimensional space of the cold region tunnel are marked to generate the distribution map of the critical area of ​​frost heave risk.

[0019] Based on the freezing front advance velocity of each segment in the distribution map of critical areas for frost heave risk, a time-series priority ranking is established by calculating the dynamic equilibrium response time between the freezing front and the unfrozen water interface, resulting in a heating response sequence with time constraints, including:

[0020] Based on the freezing front advance speed and the convergence intensity, the increase in pore ice volume caused by the freezing front advance per unit time and the volumetric flow rate of unfrozen water migration and replenishment in the convergence area are calculated. By judging the matching relationship between the pore ice volume increase and the volumetric flow rate, the critical water replenishment rate required to maintain the stable advance of the freezing front is identified as the assessment value of the continuous water supply capacity of each risk section.

[0021] Find the temperature increase required to reduce the advance speed to zero on the response characteristic curve of the freezing front advance speed as a function of temperature. Find the temperature increase required to reduce the water supply rate to the point where it can no longer sustain the advance of the freezing front on the decay characteristic curve of the water supply capacity assessment value as a function of temperature. Select the larger of the two temperature increase values ​​as the heating temperature control target. Calculate the minimum continuous heating time required to reach the heating temperature control target and obtain the dynamic equilibrium response time of each risk zone.

[0022] The remaining available time for each risk zone from the current moment to the end of the dynamic equilibrium response time is normalized and the reciprocal is taken. This reciprocal is then combined with the water supply capacity assessment value to obtain a comprehensive threat level score.

[0023] Based on the comprehensive threat level score, each risk segment is prioritized according to time sequence to generate the heating response sequence with time dimension constraints.

[0024] Based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the dominant heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracing the evolution trajectory of the temperature isosurface. The heat transfer path network topology is obtained as follows:

[0025] A set of spatial points corresponding to the same temperature value is extracted from the three-dimensional distribution of the temperature gradient to construct a temperature isosurface. The spatial position change of the temperature isosurface in a continuous time series is tracked. The spatial position change is vector decomposed in the direction perpendicular to the advance of the freezing front and in the direction parallel to the advance of the freezing front, respectively. The displacement in each direction per unit time is calculated to obtain the moving velocity component of the temperature isosurface in the direction perpendicular to the advance of the freezing front and the diffusion velocity component in the direction parallel to the advance of the freezing front.

[0026] The spatial region where the moving velocity component is greater than the diffusion velocity component is identified as the heat transfer region, and the spatial region where the diffusion velocity component is greater than the moving velocity component is identified as the heat dissipation region.

[0027] In the heat transfer region, the shortest distance paths between adjacent temperature isosurfaces are connected along the temperature gradient direction to obtain the heat transfer path. In the heat dissipation region, the lateral expansion trajectory of the temperature isosurfaces is traced along the direction perpendicular to the temperature gradient to obtain the heat dissipation path. The path direction and branch connection relationship of the heat transfer path, as well as the diffusion direction and coverage of the heat dissipation path, are extracted. Topological mapping is performed in a three-dimensional spatial coordinate system to obtain the network topology of the transfer path, which includes the spatial location of nodes, path connection method, and path type attributes.

[0028] Within the heat transfer region, the shortest distance paths between adjacent temperature isosurfaces are connected along the temperature gradient direction to obtain the heat transfer path. Within the heat dissipation region, the lateral expansion trajectory of the temperature isosurfaces is traced along the direction perpendicular to the temperature gradient to obtain the heat dissipation path. The path orientation and branch connection relationships of the heat transfer path, as well as the diffusion direction and coverage of the heat dissipation path, are extracted, including:

[0029] In the heat advantage transfer area, adjacent temperature isosurfaces are extracted layer by layer along the temperature gradient direction from the heating source to the freezing front. The Euclidean distance between each spatial point on each layer of temperature isosurface and each spatial point on the next layer of temperature isosurface is calculated. The line connecting the point pair corresponding to the minimum value of the Euclidean distance is selected as a single transfer path.

[0030] All the single-segment transmission paths are connected in descending order of temperature to form a complete heat dominance transmission path from the heating source to the freezing front. The path direction of the heat dominance transmission path in three-dimensional space and the branch connection relationship generated when the heat dominance transmission path encounters the temperature field bifurcation position are extracted.

[0031] An initial temperature isosurface is selected within the heat dissipation region. The boundary trajectory of the initial temperature isosurface as it expands outward over time is traced in a plane perpendicular to the temperature gradient direction. The spatial coordinates of the boundary trajectory at different times are recorded. The spatial coordinates of the boundary trajectory at each time are connected to form the lateral expansion trajectory of the temperature isosurface as the heat dissipation path.

[0032] The main diffusion direction vector of the lateral expansion trajectory is calculated as the diffusion direction of the heat dissipation path relative to the temperature gradient direction, and the area of ​​the spatial region enclosed by the lateral expansion trajectory is calculated as the coverage range of the heat dissipation path in three-dimensional space.

[0033] By quantifying the contribution of unit heating power to temperature rise and energy attenuation ratio when it is transferred to the freezing front along the current path, the heat utilization efficiency coefficients at different heating locations in each section are calculated, including:

[0034] The current path from the heating source to the freezing front is divided into multiple calculation segments according to the path length. A unit heating power is applied at the beginning of each calculation segment. The temperature change caused by the unit heating power at the boundary of each calculation segment during the transfer of the unit heating power to the freezing front along the current path is tracked and the ratio of the temperature change to the unit heating power is calculated as the temperature increase contribution of the unit heating power when it is transferred to the freezing front in the current calculation segment.

[0035] In each calculation segment of the heat transfer path, the difference between the input energy at the beginning of the segment and the output energy at the end of the segment is calculated, and the ratio of the difference to the input energy is used as the heat transfer energy attenuation ratio. In each calculation segment of the heat dissipation path, the energy dissipated by the unit heating power in the diffusion direction of the heat dissipation path is calculated, and the ratio of the dissipated energy to the input energy is used as the dissipation energy attenuation ratio. The heat transfer energy attenuation ratio and the dissipation energy attenuation ratio are combined as the energy attenuation ratio of the unit heating power in the current calculation segment.

[0036] By combining the temperature increase contribution with the energy attenuation ratio, the heat utilization efficiency coefficient of different heating positions in each calculation segment of the transmission path network topology is obtained.

[0037] Based on the heating response sequence and the heat utilization efficiency coefficient, a dynamic heating source spatial configuration scheme for frost heave prevention in cold-region tunnels is generated, including:

[0038] Extract the variation law of the freezing front advance velocity at each heating position under different heating power from the heating response sequence, identify the heating power critical point where the freezing front advance velocity changes from a positive value to a zero value or a negative value as the minimum antifreeze power requirement for each heating position, and associate the minimum antifreeze power requirement with the spatial coordinates of the heating position to form the spatial distribution of antifreeze power requirement;

[0039] High-efficiency heating sections with heat utilization efficiency coefficients exceeding a preset threshold are selected from the heat utilization efficiency coefficients. Within the high-efficiency heating sections, the locations of heating sources are determined according to the spatial distribution of the antifreeze power demand. The spatial gradient change rate of the minimum antifreeze power demand is extracted from the spatial distribution of the antifreeze power demand. The locations where the spatial gradient change rate is greater than a preset change rate are identified as densely distributed heating source areas. The density of heating sources is increased within the densely distributed heating source areas.

[0040] By combining the location of the heating sources in the high-efficiency heating section, the minimum antifreeze power requirement corresponding to the location of the heating sources, and the density of the heating sources in the densely distributed area, a dynamic spatial configuration scheme for antifreeze heave of heating sources in cold-region tunnels is generated, which includes the spatial location coordinates of the heating sources, the power configuration of the heating sources, and the spacing between the heating sources.

[0041] A second aspect of the present invention provides a dynamic heating system for preventing frost heave in cold-region tunnels based on multi-physics field coupling, comprising:

[0042] The first unit is used to obtain the temperature gradient, water content distribution, and freezing front advance speed at different depths in cold region tunnels.

[0043] The second unit is used to perform coupled analysis of the temperature gradient and the moisture content distribution to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and a distribution map of the critical area of ​​frost heave risk is obtained.

[0044] The third unit is used to establish a time priority sorting by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface based on the freezing front advance speed of each segment in the distribution map of the critical area of ​​freezing heave risk, and to obtain a heating response sequence with time dimension constraints.

[0045] The fourth unit is used to identify the dominant heat transfer path and heat dissipation path between the heating source and the freezing front based on the distribution law of the temperature gradient in the three-dimensional space of the cold region tunnel by tracking the evolution trajectory of the temperature isosurface, and to obtain the heat transfer path network topology.

[0046] The fifth unit is used to calculate the heat utilization efficiency coefficient of different heating positions in each section by quantifying the contribution of temperature rise and energy attenuation ratio when a unit heating power is transferred to the freezing front along the current path in each path of the transmission path network topology.

[0047] The sixth unit is used to generate a dynamic heating source spatial configuration scheme for cold-region tunnel frost heave prevention based on the heating response sequence and the heat utilization efficiency coefficient.

[0048] A third aspect of the embodiments of the present invention,

[0049] An electronic device is provided, comprising:

[0050] processor;

[0051] Memory used to store processor-executable instructions;

[0052] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0053] Fourth aspect of the present invention,

[0054] A computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0055] The beneficial effects of this application are as follows:

[0056] By coupling the analysis of temperature gradient and moisture content distribution, a potential field driving moisture migration is formed, enabling accurate identification of risk zones for continuous moisture migration and accumulation towards the freezing front. This effectively avoids the problem of inaccurate judgment of frost heave risk areas in traditional methods. Based on the freezing front advance velocity, the dynamic equilibrium response time is calculated and a time-series priority ranking is established, forming a heating response sequence with time constraints. This solves the problems of response lag and energy waste in traditional anti-frost heave heating systems. By identifying the dominant heat transfer path and heat dissipation path from the heating source to the freezing front, a heat transfer path network topology is constructed, achieving a refined description of the heat transfer process and overcoming the deficiency of incomplete consideration of heat conduction paths in traditional heating methods. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating the dynamic heating method for preventing frost heave in cold-region tunnels based on multi-physics field coupling, as described in an embodiment of the present invention.

[0058] Figure 2 This is a schematic diagram of the potential field analysis process for water migration. Detailed Implementation

[0059] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0061] Figure 1 This is a schematic flowchart of a dynamic heating method for preventing frost heave in cold-region tunnels based on multi-physics field coupling, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0062] To obtain the temperature gradient, water content distribution, and freezing front advance velocity at different depths in tunnels in cold regions;

[0063] The coupling analysis of the temperature gradient and the moisture content distribution is used to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and the distribution map of the critical area of ​​frost heave risk is obtained.

[0064] Based on the freezing front advance speed of each segment in the distribution map of critical areas of freezing heave risk, a time priority ranking is established by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface, and a heating response sequence with time dimension constraints is obtained.

[0065] Based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracking the evolution trajectory of the temperature isosurface, and the heat transfer path network topology is obtained.

[0066] In each path of the transmission path network topology, the heat utilization efficiency coefficient of different heating positions in each section is calculated by quantifying the contribution of temperature rise and energy attenuation when a unit heating power is transmitted to the freezing front along the current path.

[0067] Based on the heating response sequence and the heat utilization efficiency coefficient, a dynamic heating source spatial configuration scheme for frost heave prevention in cold region tunnels is generated.

[0068] In one optional implementation, a coupled analysis of the temperature gradient and the moisture content distribution is performed to obtain a moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, risk zones in cold-region tunnels where moisture continuously migrates and accumulates towards the freezing front are identified, resulting in a distribution map of critical areas for frost heave risk, including:

[0069] Establish a spatial discrete grid in the three-dimensional space of tunnels in cold regions;

[0070] The three-dimensional spatial components of the temperature gradient and the numerical distribution of the water content are registered and aligned on the same spatial discrete grid nodes. The water migration driving potential field is obtained by calculating the superposition effect of the adsorption potential of ice-water phase transition caused by temperature decrease and the osmotic potential generated by the water content gradient at each grid node.

[0071] The spatial gradient of potential energy value is extracted from the water migration driving potential field. The direction of the spatial gradient is taken as the preferred water migration direction, and the magnitude of the spatial gradient is taken as the water migration driving intensity, so as to obtain a potential field distribution vector field containing direction information and intensity information.

[0072] In the potential field distribution vector field, trace the potential energy gradient streamlines along the direction of potential energy reduction from the region with water content greater than the preset reference value to the region with temperature below the freezing point, and identify the convergence region where multiple potential energy gradient streamlines intersect in space.

[0073] Calculate the convergence density and convergence intensity of the potential energy gradient streamlines within the convergence region, and mark the convergence region where the convergence density exceeds a preset density threshold as a risk zone for continuous water migration and accumulation.

[0074] The location coordinates, spatial range, and convergence intensity of the risk section in the three-dimensional space of the cold region tunnel are marked to generate the distribution map of the critical area of ​​frost heave risk.

[0075] like Figure 2 As shown, the method includes:

[0076] A spatial discrete mesh is established in the three-dimensional space of the tunnel in the cold region. The tunnel is spatially discretized using a hexahedral mesh. The mesh size can be set according to the required computational accuracy. For example, the tunnel space can be divided into cubic units with a side length of 10 cm to form a complete computational domain containing the tunnel structure and the surrounding soil. For areas with complex tunnel cross-sectional shapes, local mesh refinement techniques can be used to ensure computational accuracy at the boundaries.

[0077] The temperature gradient's three-dimensional spatial components are registered and aligned with the moisture content distribution on the same spatial discrete grid nodes to obtain the temperature field data of the tunnel surrounding rock. This data can be obtained through real-time monitoring using embedded temperature sensors or through numerical simulations of thermal conduction. Similarly, moisture content distribution data is obtained through soil moisture sensors or through numerical simulations of water seepage. At each grid node, temperature and moisture content data are stored as correlated data pairs to ensure the correspondence between the two physical fields at the same spatial location.

[0078] For each grid node, the superposition effect of the adsorption potential of the ice-water phase transition caused by temperature decrease and the permeability potential generated by the water content gradient is calculated. The adsorption potential of the ice-water phase transition is proportional to the temperature difference below the freezing point and can be expressed as the temperature difference multiplied by the phase transition coefficient; the permeability potential is proportional to the water content gradient and is expressed as the water content gradient multiplied by the permeability coefficient. The components of both in each direction are added together to obtain the comprehensive water migration driving potential field. In practical calculations, the heterogeneity of different regions of the soil can be considered, and the spatial distribution of the phase transition coefficient and permeability coefficient can be corrected to improve the calculation accuracy.

[0079] For the potential energy value of each grid node, its partial derivatives in the x, y, and z directions are calculated to form a gradient vector. The direction of the gradient vector is the preferential migration direction of water, pointing towards the direction of the fastest decrease in potential energy; the magnitude of the gradient vector represents the driving force of water migration. A central difference scheme is used to calculate the spatial gradient, ensuring the stability and second-order accuracy of the calculation results. For grid points at the boundary of the computational domain, a one-sided difference method is used to calculate the gradient, ensuring the integrity of the calculation. Through this step, a potential field distribution vector field containing both directional and intensity information is obtained.

[0080] In the potential field distribution vector field, potential energy gradient streamlines are traced. Regions with a water content greater than a preset benchmark are selected as the starting points for the streamlines; for example, grid points with a water content greater than 20% can be set as starting points. Regions with temperatures below the freezing point are selected as potential endpoint regions. Starting from each starting point, the streamlines are traced step by step along the direction of decreasing potential energy, and the calculated trajectory forms the potential energy gradient streamlines. During the tracing process, the fourth-order Runge-Kutta integral method is used to solve the streamline trajectory to ensure tracing accuracy. When a streamline enters a region with a temperature below the freezing point, the endpoint position of the streamline is recorded. In this way, multiple potential energy gradient streamlines from high water content regions to low temperature regions can be obtained.

[0081] The computational domain is divided into several statistical units, and the number of streamlines passing through each unit is counted. A kernel density estimation method is used to smooth the streamline distribution, avoiding statistical errors caused by spatial discretization. For regions with dense streamlines, the statistical units are further subdivided to improve spatial resolution. This method can accurately identify spatial convergence areas of streamlines, representing locations where water accumulates.

[0082] The convergence density and convergence intensity of potential energy gradient streamlines within the convergence region are calculated. Convergence density is defined as the number of streamlines per unit volume, representing the tendency for water to converge into the region from multiple directions. Convergence intensity considers the magnitude of the potential energy gradient carried by the streamlines and is calculated as the weighted average of the potential energy gradient moduli of all streamlines passing through the region. Convergence regions where the convergence density exceeds a preset density threshold are designated as risk areas for continued water migration and accumulation. For example, a threshold of twice the average convergence density can be set to mark high-risk areas.

[0083] Information on risk sections within the three-dimensional space of tunnels in cold regions is marked to generate a distribution map of critical areas for frost heave risk. The center coordinates, spatial extent (length, width, and height), and convergence intensity value of each risk section are recorded, forming structured data. Tunnel mileage information is appended to the distribution map to facilitate accurate location of risk areas by construction and operation / maintenance personnel, enabling them to implement targeted frost heave prevention measures, such as installing insulation layers, adding drainage systems, or using antifreeze materials.

[0084] The distribution map of critical frost heave risk obtained through the above steps can provide a scientific basis for the design, construction and operation of tunnels in cold regions, enabling accurate prediction and effective prevention and control of frost heave risk, and improving the safety and durability of tunnel structures.

[0085] In one optional implementation, based on the freezing front advance velocity of each segment in the frost heave risk critical area distribution map, a time-series priority ranking is established by calculating the dynamic equilibrium response time between the freezing front and the unfrozen water interface, resulting in a heating response sequence with time-dimensional constraints, including:

[0086] Based on the freezing front advance speed and the convergence intensity, the increase in pore ice volume caused by the freezing front advance per unit time and the volumetric flow rate of unfrozen water migration and replenishment in the convergence area are calculated. By judging the matching relationship between the pore ice volume increase and the volumetric flow rate, the critical water replenishment rate required to maintain the stable advance of the freezing front is identified as the assessment value of the continuous water supply capacity of each risk section.

[0087] Find the temperature increase required to reduce the advance speed to zero on the response characteristic curve of the freezing front advance speed as a function of temperature. Find the temperature increase required to reduce the water supply rate to the point where it can no longer sustain the advance of the freezing front on the decay characteristic curve of the water supply capacity assessment value as a function of temperature. Select the larger of the two temperature increase values ​​as the heating temperature control target. Calculate the minimum continuous heating time required to reach the heating temperature control target and obtain the dynamic equilibrium response time of each risk zone.

[0088] The remaining available time for each risk zone from the current moment to the end of the dynamic equilibrium response time is normalized and the reciprocal is taken. This reciprocal is then combined with the water supply capacity assessment value to obtain a comprehensive threat level score.

[0089] Based on the comprehensive threat level score, each risk segment is prioritized according to time sequence to generate the heating response sequence with time dimension constraints.

[0090] Based on the advance velocity of the freezing front in each segment of the critical zone distribution map for frost heave risk, the increase in pore ice volume caused by the advance of the freezing front per unit time can be calculated. Simultaneously, the matching relationship between these two parameters is determined by measuring the volumetric flow rate of unfrozen water migration and replenishment in the convergence area. When the increase in pore ice volume generated by the advance of the freezing front equals the volumetric flow rate of unfrozen water migration and replenishment, the freezing front will maintain stable advance. The water replenishment rate under this critical state is defined as the assessment value of the sustainable water supply capacity of this risk segment.

[0091] Specifically, the following calculation method can be used: For each segment within the frost heave risk area, soil porosity, ice density, and the advance velocity of the freezing front are measured to calculate the newly added ice volume per unit time. Simultaneously, the volumetric flow rate of water recharge in the convergence area is calculated by measuring the flow rate and flow area of ​​unfrozen water. When these two volumetric changes reach equilibrium, the freezing front will maintain a stable advance.

[0092] This study analyzes the response characteristics of the freezing front's advance velocity to temperature changes, establishes a curve showing the relationship between temperature and the freezing front's advance velocity through experiments, and identifies the temperature increase required to reduce the advance velocity to zero. This temperature increase marks the critical point at which the freezing process will be completely suppressed. Simultaneously, the continuous water supply capacity also decreases with increasing temperature; therefore, a decay characteristic curve of the water supply rate as a function of temperature needs to be established. The critical temperature increase required to reduce the water supply rate to a level that makes it impossible to sustain the freezing front's advance is then identified on this curve.

[0093] Comparing the two temperature increases mentioned above, the larger value is selected as the heating temperature control target. This ensures that both the freezing process is suppressed and the rate of moisture migration is controlled. Based on the selected heating temperature control target, the minimum continuous heating time required to reach this temperature is calculated using a heat conduction model, which is the dynamic equilibrium response time for this risk zone.

[0094] In practical applications, temperature sensor arrays can be used to monitor temperature changes on the freezing front, and humidity sensors buried at different depths can be used to monitor moisture migration. The heater power and heating time can be set based on the calculated temperature rise and minimum continuous heating duration.

[0095] To allocate heating resources rationally, priority needs to be assigned to each risk zone. The remaining available time for each zone from the current moment until the end of the dynamic equilibrium response time is calculated, normalized, and the reciprocal is taken to obtain a time urgency index. This index is then combined with the water sustainability assessment value to obtain a comprehensive threat level score.

[0096] For example, assuming a section has a short remaining available time but a strong water supply capacity, its overall threat level score will be high, requiring priority for heating treatment. A weighted average method can be used to combine the time urgency index and the water supply capacity assessment value: Overall Threat Level Score = α × Normalized Time Urgency Index + β × Normalized Water Supply Capacity Assessment Value. Here, α and β are weighting coefficients that can be adjusted according to actual engineering needs.

[0097] Based on the comprehensive threat level score of each risk zone, a time-series priority ranking is performed to generate a heating response sequence with time constraints. This sequence guides heating equipment to precisely heat each risk zone according to priority, ensuring that the overall risk of frost heave is minimized under limited resource conditions.

[0098] In one optional implementation, based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the dominant heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracing the evolution trajectory of the temperature isosurface, resulting in the heat transfer path network topology including:

[0099] A set of spatial points corresponding to the same temperature value is extracted from the three-dimensional distribution of the temperature gradient to construct a temperature isosurface. The spatial position change of the temperature isosurface in a continuous time series is tracked. The spatial position change is vector decomposed in the direction perpendicular to the advance of the freezing front and in the direction parallel to the advance of the freezing front, respectively. The displacement in each direction per unit time is calculated to obtain the moving velocity component of the temperature isosurface in the direction perpendicular to the advance of the freezing front and the diffusion velocity component in the direction parallel to the advance of the freezing front.

[0100] The spatial region where the moving velocity component is greater than the diffusion velocity component is identified as the heat transfer region, and the spatial region where the diffusion velocity component is greater than the moving velocity component is identified as the heat dissipation region.

[0101] In the heat transfer region, the shortest distance paths between adjacent temperature isosurfaces are connected along the temperature gradient direction to obtain the heat transfer path. In the heat dissipation region, the lateral expansion trajectory of the temperature isosurfaces is traced along the direction perpendicular to the temperature gradient to obtain the heat dissipation path. The path direction and branch connection relationship of the heat transfer path, as well as the diffusion direction and coverage of the heat dissipation path, are extracted. Topological mapping is performed in a three-dimensional spatial coordinate system to obtain the network topology of the transfer path, which includes the spatial location of nodes, path connection method, and path type attributes.

[0102] A temperature monitoring network is established in the three-dimensional space of a tunnel in a cold region. Multiple temperature sensors are deployed to collect temperature data at different locations and times. Based on the collected temperature data, a temperature field distribution model of the entire tunnel space is constructed using three-dimensional interpolation methods. Kriging interpolation or radial basis function interpolation can effectively process discrete point temperature data and generate a continuous temperature field model.

[0103] After obtaining the temperature field distribution, extract the set of spatial points corresponding to the same temperature value, construct a temperature isosurface, and select a series of temperature thresholds T1, T2...T n For each threshold T i Find all pairs of elements in three-dimensional space that satisfy T(x, y, z) = T i The points (x, y, z) form a closed surface, namely the temperature isosurface. The Marching Cubes algorithm can efficiently extract these isosurfaces from the three-dimensional temperature field and represent them as a triangular mesh.

[0104] To track the temporal evolution of temperature isosurfaces, in a continuous time series t1, t2...t m The positional changes of each isosurface are recorded. For time point tj, the temperature isosurface S(T) is... i , t j ) and time point t j+1 The corresponding temperature isosurface S(T) i , t j+1 Calculate the spatial displacement between two isosurfaces. Select a set of feature points P1, P2...P on the isosurfaces. k The displacement vector of each feature point within the time interval is calculated using the nearest point matching algorithm.

[0105] The displacement vector needs to be decomposed into a component perpendicular to the direction of the freezing front and a component parallel to the direction of the freezing front. Let the normal vector of the freezing front be n. For feature point P... i The displacement vector d has a vertical component d. v =(d·n)n, where the parallel component is d. p =dd vCalculate the displacement per unit time, with the vertical velocity component being v. v =|d v | / (t j +1-t j The diffusion velocity component in the parallel direction is v. p =|d p | / (t j +1-t j ).

[0106] Based on the calculated velocity components, the regions of dominant heat transfer and heat dissipation are identified. When the moving velocity component is greater than the diffusion velocity component (v... v >v p When the diffusion velocity component is greater than the moving velocity component (v), the region is identified as a region with a dominant heat transfer; conversely, when the diffusion velocity component is greater than the moving velocity component (v), the region is identified as a region with a dominant heat transfer. p >v v When this occurs, the region is identified as a heat dissipation region.

[0107] Within the region of advantageous heat transfer, the heat transfer path is traced along the temperature gradient direction. The temperature gradient direction can be obtained by calculating the gradient vector ▽T(x, y, z) of the temperature field. Starting from the heating source, the shortest path is found between adjacent isosurfaces along the gradient direction. In practice, the A* search algorithm is used, with the temperature gradient direction as the heuristic function, to search for the optimal path between adjacent isosurfaces and then concatenate these paths to form a complete advantageous heat transfer path.

[0108] Within the heat dissipation region, the lateral expansion trajectory of the temperature isosurface needs to be tracked. Perpendicular to the temperature gradient direction, curved slices are made along the isosurface, and the changes in the contour line of the isosurface on the slice plane over time are recorded. By comparing the changes in the area and shape of the contour line at different times, the lateral diffusion direction and rate of heat are determined, thus obtaining the heat dissipation path.

[0109] Heat transfer and dissipation paths are integrated into a three-dimensional spatial coordinate system to construct a network topology for the transfer paths. For each path, its spatial coordinates, path type attribute (superior transfer or dissipation), and connection relationships with other paths are recorded. Connection points between paths are defined as network nodes, and the spatial location of each node and the number of paths it connects to are recorded.

[0110] In an application example, during a tunnel project in a cold region, the heat transfer network identified using the aforementioned method revealed several dominant heat transfer paths within the tunnel's surrounding rock. These paths often follow rock fissures or aquifers. Simultaneously, multiple heat dissipation zones were identified at the tunnel's top and sidewalls, areas requiring focused attention for frost protection measures. Based on this heat transfer network topology, the tunnel insulation design was optimized, resulting in more targeted heat source placement, improved energy efficiency, and effective prevention of frost damage.

[0111] This method of identifying heat transfer path networks can provide a scientific basis for the thermal design and frost damage prevention of tunnels in cold regions, and has important engineering application value.

[0112] In one optional implementation, a heat-dominant transfer path is obtained by connecting the shortest distance paths between adjacent temperature isosurfaces along the temperature gradient direction within the heat-dominant transfer region. A heat-dissipation path is obtained by tracing the lateral expansion trajectory of the temperature isosurfaces along a direction perpendicular to the temperature gradient within the heat-dissipation region. Extracting the path orientation and branch connections of the heat-dominant transfer path, as well as the diffusion direction and coverage of the heat-dissipation path, includes:

[0113] In the heat advantage transfer area, adjacent temperature isosurfaces are extracted layer by layer along the temperature gradient direction from the heating source to the freezing front. The Euclidean distance between each spatial point on each layer of temperature isosurface and each spatial point on the next layer of temperature isosurface is calculated. The line connecting the point pair corresponding to the minimum value of the Euclidean distance is selected as a single transfer path.

[0114] All the single-segment transmission paths are connected in descending order of temperature to form a complete heat dominance transmission path from the heating source to the freezing front. The path direction of the heat dominance transmission path in three-dimensional space and the branch connection relationship generated when the heat dominance transmission path encounters the temperature field bifurcation position are extracted.

[0115] An initial temperature isosurface is selected within the heat dissipation region. The boundary trajectory of the initial temperature isosurface as it expands outward over time is traced in a plane perpendicular to the temperature gradient direction. The spatial coordinates of the boundary trajectory at different times are recorded. The spatial coordinates of the boundary trajectory at each time are connected to form the lateral expansion trajectory of the temperature isosurface as the heat dissipation path.

[0116] The main diffusion direction vector of the lateral expansion trajectory is calculated as the diffusion direction of the heat dissipation path relative to the temperature gradient direction, and the area of ​​the spatial region enclosed by the lateral expansion trajectory is calculated as the coverage range of the heat dissipation path in three-dimensional space.

[0117] Temperature field data is collected from the heat flow field by deploying an array of temperature sensors in the heat transfer space to obtain three-dimensional spatial temperature distribution data. Based on the collected temperature distribution data, regions with dominant heat transfer and heat dissipation can be distinguished. The dominant heat transfer region is characterized by a significant temperature gradient, while the heat dissipation region is characterized by heat diffusing into the surrounding environment.

[0118] Within a region of dominant heat transfer, the process of extracting the dominant heat transfer path along the temperature gradient is as follows: Starting from the location of the heating source, temperature isosurfaces are extracted layer by layer in descending order of temperature. For example, when the heating source temperature is 80℃ and the freezing front temperature is 0℃, a temperature isosurface can be extracted every 5℃, resulting in a series of temperature isosurfaces such as 80℃, 75℃, 70℃...5℃, 0℃, etc.

[0119] For each pair of adjacent temperature isosurfaces, such as a 75℃ temperature isosurface and a 70℃ temperature isosurface, calculate the Euclidean distance from each point on the 75℃ temperature isosurface to all points on the 70℃ temperature isosurface. Assuming there is a point P1(x1, y1, z1) on the 75℃ temperature isosurface and a point P2(x2, y2, z2) on the 70℃ temperature isosurface, the Euclidean distance between the two points is calculated as the spatial distance between their coordinates. By comparing the distances of all point pairs, select the point pair with the minimum distance, and use the line connecting these two points as a single-segment transmission path.

[0120] Repeat the above process, starting from the heat source, calculating the single-segment heat transfer paths between each temperature isosurface layer by layer, and then connecting these single-segment paths in series according to the temperature from high to low to form a complete dominant heat transfer path from the heat source to the freezing front. For example, connecting the paths from 80℃ to 75℃, 75℃ to 70℃... 5℃ to 0℃ in series yields the complete dominant heat transfer path from the heat source to the freezing front.

[0121] In the process of extracting the dominant heat transfer path, temperature field bifurcation may occur, meaning that multiple shortest paths with similar Euclidean distances appear between two isosurfaces. In this case, the positions and connections of these bifurcation points are recorded to construct the topology of the dominant heat transfer path. For example, when heat transfer occurs from point P on the 50℃ isosurface to the 45℃ isosurface, there will be two or more arrival points Q1 and Q2, thus forming a branch of the path.

[0122] Within the heat dissipation region, the process of extracting the heat dissipation path is as follows: An initial temperature isosurface is selected, for example, a 30°C temperature isosurface as the interface where significant heat dissipation begins. In a plane perpendicular to the temperature gradient direction, the lateral expansion of this temperature isosurface over time is tracked. Specifically, at a series of times t1, t2, t3, etc., the spatial coordinates of the temperature isosurface boundary are recorded. Connecting these coordinate points forms the lateral expansion trajectory of the temperature isosurface, i.e., the heat dissipation path.

[0123] For the obtained heat dissipation path, calculate its main diffusion direction vector and the centroid coordinates of the heat dissipation path at different times. Use these centroid coordinates to fit a main diffusion direction vector, which represents the main diffusion direction of heat relative to the temperature gradient direction. For example, if the heat mainly diffuses in the positive x-axis direction, then the main diffusion direction vector will point in the positive x-axis direction.

[0124] The extent of heat dissipation paths in three-dimensional space is calculated. For each moment, the area of ​​the spatial region enclosed by the temperature isosurface boundary trajectory is calculated. This area gradually increases over time, indicating that the range of heat diffusion into the surrounding environment is expanding. By recording the coverage area at different moments, the relationship between the spatial coverage of heat dissipation and time can be obtained.

[0125] The heat transfer and dissipation paths extracted using the methods described above can visually represent the heat transfer and diffusion characteristics in three-dimensional space. For example, in geothermal harvesting systems, the optimal heat transfer path from the geothermal source to the harvesting equipment, as well as the main directions and extent of heat loss to the surrounding environment, can be identified, providing a basis for the optimized design of geothermal harvesting systems.

[0126] In practical applications, the density of temperature field data acquisition affects the accuracy of heat transfer and dissipation path extraction. Higher acquisition density provides more refined temperature gradient information, resulting in more accurate heat transfer and dissipation paths. However, excessively high acquisition density also increases computational load, requiring a balance between accuracy and efficiency. Generally, the density of temperature acquisition points can be adjusted according to the needs of the actual application, with denser placement of points in key areas to improve the accuracy of heat path extraction.

[0127] In one optional implementation, the heat utilization efficiency coefficients for different heating locations in each segment are calculated by quantifying the contribution of temperature rise and energy attenuation when a unit heating power is transferred along the current path to the freezing front, including:

[0128] The current path from the heating source to the freezing front is divided into multiple calculation segments according to the path length. A unit heating power is applied at the beginning of each calculation segment. The temperature change caused by the unit heating power at the boundary of each calculation segment during the transfer of the unit heating power to the freezing front along the current path is tracked and the ratio of the temperature change to the unit heating power is calculated as the temperature increase contribution of the unit heating power when it is transferred to the freezing front in the current calculation segment.

[0129] In each calculation segment of the heat transfer path, the difference between the input energy at the beginning of the segment and the output energy at the end of the segment is calculated, and the ratio of the difference to the input energy is used as the heat transfer energy attenuation ratio. In each calculation segment of the heat dissipation path, the energy dissipated by the unit heating power in the diffusion direction of the heat dissipation path is calculated, and the ratio of the dissipated energy to the input energy is used as the dissipation energy attenuation ratio. The heat transfer energy attenuation ratio and the dissipation energy attenuation ratio are combined as the energy attenuation ratio of the unit heating power in the current calculation segment.

[0130] By combining the temperature increase contribution with the energy attenuation ratio, the heat utilization efficiency coefficient of different heating positions in each calculation segment of the transmission path network topology is obtained.

[0131] The calculation area is divided according to the geological model. The heat transfer path from the heating source to the freezing front is divided into several calculation sections according to the path length. For example, if the total length of the transfer path is 100 meters, it can be divided into 10 sections, each 10 meters long. Each section is defined by a starting boundary and a ending boundary, forming a discretized calculation basis.

[0132] When calculating the contribution of heating to temperature rise, a unit heating power, such as 1000 watts, is applied at the beginning of each calculation segment. The propagation of this power along the transmission path to the freezing front is tracked using heat conduction simulation. Temperature changes are recorded at the boundaries of each segment. For example, if the temperature rises by 2.5°C at the end boundary of the first segment, the contribution of the unit heating power to temperature rise in that segment is 0.0025°C / watt. This method can visually reflect the effect of heating at different locations on the temperature of the target freezing front.

[0133] The energy attenuation ratio calculation is divided into two parts: dominant transfer energy attenuation and dissipated energy attenuation. For each calculation segment of the dominant transfer path, the input energy per unit heating power at the start of the segment and the output energy at the end of the segment are recorded. For example, in a certain segment, the input energy is 800 joules, and the output energy at the end is 720 joules, then the difference is 80 joules, and the dominant transfer energy attenuation ratio is 80 / 800 = 10%. This indicates that 10% of the energy in this segment failed to transfer along the dominant path.

[0134] For each calculation segment of the heat dissipation path, the energy dissipated per unit heating power in the direction of diffusion to the surrounding environment is calculated. Using a heat diffusion model, the energy value diffused outward from the dominant transfer path is determined. For example, if in a certain segment, the energy dissipated to the surrounding environment is 150 joules, and the input energy is 800 joules, then the energy dissipation attenuation ratio is 150 / 800 = 18.75%. This reflects the degree of heat loss in non-target directions during the transfer process.

[0135] By combining the energy attenuation ratio of the advantageous transfer and the energy attenuation ratio of the dissipated energy, we obtain the total energy attenuation ratio of the unit heating power in the current calculation segment. In the example above, the total energy attenuation ratio is 10% + 18.75% = 28.75%, indicating that 28.75% of the energy in this segment was not effectively transferred to the next segment.

[0136] Based on the contribution of temperature increase and the proportion of energy attenuation, the heat utilization efficiency coefficient for different heating positions in each calculation section can be obtained. The combination method can use a weighted average, where the contribution of temperature increase has a weight of 0.6 and the proportion of energy attenuation has a weight of 0.4. For the aforementioned example, the heat utilization efficiency coefficient can be expressed as 0.6 × 0.0025 - (0.4 × 0.2875) = 0.0015 - 0.115 = -0.1135.

[0137] In practical applications, simulations were conducted for a geothermal field, assuming the heat transfer path was divided into five segments. The first segment is closest to the heat source, and the fifth segment is adjacent to the freezing front. Experimental data showed that when a unit power was applied at the beginning of the first segment, only about 20% of the heat was transferred to the freezing front; while when the same power was applied at the beginning of the fourth segment, about 75% of the heat reached the freezing front. By systematically calculating the heat utilization efficiency coefficient of each segment, it can be concluded that the segment closer to the freezing front has a higher heat utilization efficiency.

[0138] This method can also take into account geological heterogeneity. Heat conduction characteristics vary significantly across different geological structures. By adjusting the calculation parameters for each section, it can adapt to the needs of thermal efficiency assessment under complex geological conditions. For example, in fault zones or high-permeability areas, the heat transfer rate is faster and the energy attenuation ratio is lower; while in low-permeability areas, the weight of the energy attenuation ratio needs to be increased for the corresponding sections.

[0139] In summary, by quantifying the contribution of unit heating power to temperature increase and the proportion of energy attenuation along the transmission path, the heat utilization efficiency coefficient of different heating locations in each section can be scientifically calculated, providing a theoretical basis and technical support for optimizing heating locations in geothermal resource development.

[0140] In one optional implementation, generating a dynamic heating source spatial configuration scheme for frost heave prevention in cold-region tunnels, based on the heating response sequence and the heat utilization efficiency coefficient, includes:

[0141] Extract the variation law of the freezing front advance velocity at each heating position under different heating power from the heating response sequence, identify the heating power critical point where the freezing front advance velocity changes from a positive value to a zero value or a negative value as the minimum antifreeze power requirement for each heating position, and associate the minimum antifreeze power requirement with the spatial coordinates of the heating position to form the spatial distribution of antifreeze power requirement;

[0142] High-efficiency heating sections with heat utilization efficiency coefficients exceeding a preset threshold are selected from the heat utilization efficiency coefficients. Within the high-efficiency heating sections, the locations of heating sources are determined according to the spatial distribution of the antifreeze power demand. The spatial gradient change rate of the minimum antifreeze power demand is extracted from the spatial distribution of the antifreeze power demand. The locations where the spatial gradient change rate is greater than a preset change rate are identified as densely distributed heating source areas. The density of heating sources is increased within the densely distributed heating source areas.

[0143] By combining the location of the heating sources in the high-efficiency heating section, the minimum antifreeze power requirement corresponding to the location of the heating sources, and the density of the heating sources in the densely distributed area, a dynamic spatial configuration scheme for antifreeze heave of heating sources in cold-region tunnels is generated, which includes the spatial location coordinates of the heating sources, the power configuration of the heating sources, and the spacing between the heating sources.

[0144] The heating response sequence was analyzed to extract the variation pattern of the freezing front advance velocity at each heating location under different heating powers. Specifically, time-series analysis was performed on the temperature data of each monitoring point to calculate the position of the freezing front (i.e., the zero-degree isotherm) at each time point. The freezing front advance velocity was obtained by dividing the difference in freezing front position between adjacent time points by the time interval. For example, if at heating location P1, the freezing front velocity is -2 mm / h (negative values ​​indicate freezing front retreat) at 100W, -0.5 mm / h at 80W, 0.2 mm / h at 60W, and 1 mm / h at 40W, then through interpolation, the critical power for the freezing front advance velocity at this location to change from positive to zero or negative is approximately 65W. This value is the minimum antifreeze power requirement for this location.

[0145] Repeat the above calculations for multiple heating locations within the space surrounding the tunnel to obtain a series of spatial points and their corresponding minimum antifreeze power requirements, constituting a spatial distribution of antifreeze power requirements. For example, a dataset like {(x1, y1, z1, P1), (x2, y2, z2, P2)...(xn, yn, zn, Pn)} can be obtained, where (xi, yi, zi) represent spatial coordinates and Pi represents the minimum antifreeze power requirement at that location.

[0146] The distribution of heat utilization efficiency coefficients (HQCs) is analyzed. HQCs represent the antifreeze effect produced per unit of heating power and are typically related to factors such as soil properties, moisture content, and burial depth. Based on pre-acquired HQC data, a reasonable threshold (e.g., 0.75) is set, and areas with HQCs exceeding this threshold are selected as high-efficiency heating sections. For example, the HQC of the area 0.5 to 1.5 meters below the tunnel arch is 0.82, the HQC of the area within 1 meter of the tunnel side walls is 0.78, while the HQC of the area below the tunnel floor is only 0.65. Therefore, the area near the arch and side walls is identified as a high-efficiency heating section.

[0147] Within a defined high-efficiency heating zone, heating sources are arranged according to the spatial distribution of antifreeze power demand, and the spatial gradient rate of change of antifreeze power demand is calculated. The spatial gradient rate of change is calculated as follows: In three-dimensional space, arbitrarily select a direction (e.g., the x-direction), calculate the power demand difference between adjacent measuring points, and then calculate the rate of change of this difference along space. If the spatial gradient rate of change in a certain region exceeds a preset threshold (e.g., 10 W / m²), the calculation is performed. 2 If the area is identified as a region with a dense distribution of heating sources, then the area will be identified as such.

[0148] Taking the tunnel arch as an example, if a measuring point is taken every 0.5 meters along the tunnel's longitudinal direction, and the measured power demands are 80W, 85W, 95W, 120W, and 130W respectively, then the power demand gradient between the third and fourth points is (120-95) / 0.5 = 50W / m, which is significantly higher than the preset threshold. Therefore, this area is marked as a densely distributed heating source area. In the standard area, the spacing between heating sources can be set to 1.5 meters, while in the densely distributed area, the spacing can be shortened to 0.8 meters to cope with the larger power demand changes in this area.

[0149] The system integrates the determined locations of heating sources within the high-efficiency heating zone, the corresponding minimum anti-freeze power requirements, and the spacing between heating sources in densely populated areas to generate a complete spatial configuration scheme. This scheme specifically includes: the three-dimensional spatial coordinates of the heating sources (x, y, z), the power value configured for each heating source (with a margin of 10% to 20% based on the minimum anti-freeze power requirement at that location), and the spacing between heating sources in different areas.

[0150] In practical applications, such as a tunnel project on the Qinghai-Tibet Plateau, the configuration scheme obtained by using the above method shows that: a heating source is placed every 1.5 meters along the longitudinal direction of the tunnel at a distance of 0.8 meters below the tunnel arch, with a power configuration of 85W; in the area with deeper permafrost at the tunnel entrance and exit sections (longitudinal mileage K0+020 to K0+050), the spacing between heating sources is reduced to 0.7 meters, and the power is increased to 120W; in the area of ​​the tunnel side walls, the heating sources are placed 1 meter away from the tunnel wall, with a spacing of 2 meters, and a power configuration of 75W.

[0151] This configuration scheme fully considers the distribution characteristics and thermal effects of permafrost, and achieves a reasonable allocation of heating resources. Through field verification, this scheme saves about 30% of the total power consumption compared with the traditional uniform layout method. At the same time, it effectively prevents tunnel structure deformation caused by frost heave and ensures the safe and stable operation of the tunnel structure in winter.

[0152] Furthermore, in practical applications, the heating scheme can be dynamically adjusted according to the tunnel construction progress and seasonal changes. For example, in the early stages of winter, some heating sources can be activated initially, and as the temperature further decreases, the power can be gradually increased or more heating sources can be activated to achieve the optimal balance between energy consumption and antifreeze effect.

[0153] This invention relates to a dynamic heating system for preventing frost heave in cold-region tunnels based on multi-physics field coupling. The system includes:

[0154] The first unit is used to obtain the temperature gradient, water content distribution, and freezing front advance speed at different depths in cold region tunnels.

[0155] The second unit is used to perform coupled analysis of the temperature gradient and the moisture content distribution to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and a distribution map of the critical area of ​​frost heave risk is obtained.

[0156] The third unit is used to establish a time priority sorting by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface based on the freezing front advance speed of each segment in the distribution map of the critical area of ​​freezing heave risk, and to obtain a heating response sequence with time dimension constraints.

[0157] The fourth unit is used to identify the dominant heat transfer path and heat dissipation path between the heating source and the freezing front based on the distribution law of the temperature gradient in the three-dimensional space of the cold region tunnel by tracking the evolution trajectory of the temperature isosurface, and to obtain the heat transfer path network topology.

[0158] The fifth unit is used to calculate the heat utilization efficiency coefficient of different heating positions in each section by quantifying the contribution of temperature rise and energy attenuation ratio when a unit heating power is transferred to the freezing front along the current path in each path of the transmission path network topology.

[0159] The sixth unit is used to generate a dynamic heating source spatial configuration scheme for cold-region tunnel frost heave prevention based on the heating response sequence and the heat utilization efficiency coefficient.

[0160] A third aspect of the present invention provides an electronic device, comprising:

[0161] processor;

[0162] Memory used to store processor-executable instructions;

[0163] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0164] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0165] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0166] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A dynamic heating method for preventing frost heave in cold-region tunnels based on multi-physics field coupling, characterized in that, include: To obtain the temperature gradient, water content distribution, and freezing front advance velocity at different depths in tunnels in cold regions; The coupling analysis of the temperature gradient and the moisture content distribution is used to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and the distribution map of the critical area of ​​frost heave risk is obtained. Based on the freezing front advance speed of each segment in the distribution map of critical areas of freezing heave risk, a time priority ranking is established by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface, and a heating response sequence with time dimension constraints is obtained. Based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracking the evolution trajectory of the temperature isosurface, and the heat transfer path network topology is obtained. In each path of the transmission path network topology, the heat utilization efficiency coefficient of different heating positions in each section is calculated by quantifying the contribution of temperature rise and energy attenuation when a unit heating power is transmitted to the freezing front along the current path. Based on the heating response sequence and the heat utilization efficiency coefficient, a dynamic heating source spatial configuration scheme for frost heave prevention in cold region tunnels is generated.

2. The method according to claim 1, characterized in that, The coupling analysis of the temperature gradient and the moisture content distribution yields the driving potential field for moisture migration. Based on the spatial distribution characteristics of this driving potential field, risk zones in cold-region tunnels where moisture continuously migrates and accumulates towards the freezing front are identified, resulting in a distribution map of critical areas for frost heave risk, including: Establish a spatial discrete grid in the three-dimensional space of tunnels in cold regions; The three-dimensional spatial components of the temperature gradient and the numerical distribution of the water content are registered and aligned on the same spatial discrete grid nodes. The water migration driving potential field is obtained by calculating the superposition effect of the adsorption potential of ice-water phase transition caused by temperature decrease and the osmotic potential generated by the water content gradient at each grid node. The spatial gradient of potential energy value is extracted from the water migration driving potential field. The direction of the spatial gradient is taken as the preferred water migration direction, and the magnitude of the spatial gradient is taken as the water migration driving intensity, so as to obtain a potential field distribution vector field containing direction information and intensity information. In the potential field distribution vector field, trace the potential energy gradient streamlines along the direction of potential energy reduction from the region with water content greater than the preset reference value to the region with temperature below the freezing point, and identify the convergence region where multiple potential energy gradient streamlines intersect in space. Calculate the convergence density and convergence intensity of the potential energy gradient streamlines within the convergence region, and mark the convergence region where the convergence density exceeds a preset density threshold as a risk zone for continuous water migration and accumulation. The location coordinates, spatial range, and convergence intensity of the risk section in the three-dimensional space of the cold region tunnel are marked to generate the distribution map of the critical area of ​​frost heave risk.

3. The method according to claim 2, characterized in that, Based on the freezing front advance velocity of each segment in the distribution map of critical areas for frost heave risk, a time-series priority ranking is established by calculating the dynamic equilibrium response time between the freezing front and the unfrozen water interface, resulting in a heating response sequence with time constraints, including: Based on the freezing front advance speed and the convergence intensity, the increase in pore ice volume caused by the freezing front advance per unit time and the volumetric flow rate of unfrozen water migration and replenishment in the convergence area are calculated. By judging the matching relationship between the pore ice volume increase and the volumetric flow rate, the critical water replenishment rate required to maintain the stable advance of the freezing front is identified as the assessment value of the continuous water supply capacity of each risk section. Find the temperature increase required to reduce the advance speed to zero on the response characteristic curve of the freezing front advance speed as a function of temperature. Find the temperature increase required to reduce the water supply rate to the point where it can no longer sustain the advance of the freezing front on the decay characteristic curve of the water supply capacity assessment value as a function of temperature. Select the larger of the two temperature increase values ​​as the heating temperature control target. Calculate the minimum continuous heating time required to reach the heating temperature control target and obtain the dynamic equilibrium response time of each risk zone. The remaining available time for each risk zone from the current moment to the end of the dynamic equilibrium response time is normalized and the reciprocal is taken. This reciprocal is then combined with the water supply capacity assessment value to obtain a comprehensive threat level score. Based on the comprehensive threat level score, each risk segment is prioritized according to time sequence to generate the heating response sequence with time dimension constraints.

4. The method according to claim 1, characterized in that, Based on the distribution pattern of the temperature gradient in the three-dimensional space of the cold region tunnel, the dominant heat transfer path and heat dissipation path between the heating source and the freezing front are identified by tracing the evolution trajectory of the temperature isosurface. The heat transfer path network topology is obtained as follows: A set of spatial points corresponding to the same temperature value is extracted from the three-dimensional distribution of the temperature gradient to construct a temperature isosurface. The spatial position change of the temperature isosurface in a continuous time series is tracked. The spatial position change is vector decomposed in the direction perpendicular to the advance of the freezing front and in the direction parallel to the advance of the freezing front, respectively. The displacement in each direction per unit time is calculated to obtain the moving velocity component of the temperature isosurface in the direction perpendicular to the advance of the freezing front and the diffusion velocity component in the direction parallel to the advance of the freezing front. The spatial region where the moving velocity component is greater than the diffusion velocity component is identified as the heat transfer region, and the spatial region where the diffusion velocity component is greater than the moving velocity component is identified as the heat dissipation region. In the heat transfer region, the shortest distance paths between adjacent temperature isosurfaces are connected along the temperature gradient direction to obtain the heat transfer path. In the heat dissipation region, the lateral expansion trajectory of the temperature isosurfaces is traced along the direction perpendicular to the temperature gradient to obtain the heat dissipation path. The path direction and branch connection relationship of the heat transfer path, as well as the diffusion direction and coverage of the heat dissipation path, are extracted. Topological mapping is performed in a three-dimensional spatial coordinate system to obtain the network topology of the transfer path, which includes the spatial location of nodes, path connection method, and path type attributes.

5. The method according to claim 4, characterized in that, Within the heat transfer region, the shortest distance paths between adjacent temperature isosurfaces are connected along the temperature gradient direction to obtain the heat transfer path. Within the heat dissipation region, the lateral expansion trajectory of the temperature isosurfaces is traced along the direction perpendicular to the temperature gradient to obtain the heat dissipation path. The path orientation and branch connection relationships of the heat transfer path, as well as the diffusion direction and coverage of the heat dissipation path, are extracted, including: In the heat advantage transfer area, adjacent temperature isosurfaces are extracted layer by layer along the temperature gradient direction from the heating source to the freezing front. The Euclidean distance between each spatial point on each layer of temperature isosurface and each spatial point on the next layer of temperature isosurface is calculated. The line connecting the point pair corresponding to the minimum value of the Euclidean distance is selected as a single transfer path. All the single-segment transmission paths are connected in descending order of temperature to form a complete heat dominance transmission path from the heating source to the freezing front. The path direction of the heat dominance transmission path in three-dimensional space and the branch connection relationship generated when the heat dominance transmission path encounters the temperature field bifurcation position are extracted. An initial temperature isosurface is selected within the heat dissipation region. The boundary trajectory of the initial temperature isosurface as it expands outward over time is traced in a plane perpendicular to the temperature gradient direction. The spatial coordinates of the boundary trajectory at different times are recorded. The spatial coordinates of the boundary trajectory at each time are connected to form the lateral expansion trajectory of the temperature isosurface as the heat dissipation path. The main diffusion direction vector of the lateral expansion trajectory is calculated as the diffusion direction of the heat dissipation path relative to the temperature gradient direction, and the area of ​​the spatial region enclosed by the lateral expansion trajectory is calculated as the coverage range of the heat dissipation path in three-dimensional space.

6. The method according to claim 1, characterized in that, By quantifying the contribution of unit heating power to temperature rise and energy attenuation ratio when it is transferred to the freezing front along the current path, the heat utilization efficiency coefficients at different heating locations in each section are calculated, including: The current path from the heating source to the freezing front is divided into multiple calculation segments according to the path length. A unit heating power is applied at the beginning of each calculation segment. The temperature change caused by the unit heating power at the boundary of each calculation segment during the transfer of the unit heating power to the freezing front along the current path is tracked and the ratio of the temperature change to the unit heating power is calculated as the temperature increase contribution of the unit heating power when it is transferred to the freezing front in the current calculation segment. In each calculation segment of the heat transfer path, the difference between the input energy at the beginning of the segment and the output energy at the end of the segment is calculated, and the ratio of the difference to the input energy is used as the heat transfer energy attenuation ratio. In each calculation segment of the heat dissipation path, the energy dissipated by the unit heating power in the diffusion direction of the heat dissipation path is calculated, and the ratio of the dissipated energy to the input energy is used as the dissipation energy attenuation ratio. The heat transfer energy attenuation ratio and the dissipation energy attenuation ratio are combined as the energy attenuation ratio of the unit heating power in the current calculation segment. By combining the temperature increase contribution with the energy attenuation ratio, the heat utilization efficiency coefficient of different heating positions in each calculation segment of the transmission path network topology is obtained.

7. The method according to claim 1, characterized in that, Based on the heating response sequence and the heat utilization efficiency coefficient, a dynamic heating source spatial configuration scheme for frost heave prevention in cold-region tunnels is generated, including: Extract the variation law of the freezing front advance velocity at each heating position under different heating power from the heating response sequence, identify the heating power critical point where the freezing front advance velocity changes from a positive value to a zero value or a negative value as the minimum antifreeze power requirement for each heating position, and associate the minimum antifreeze power requirement with the spatial coordinates of the heating position to form the spatial distribution of antifreeze power requirement; High-efficiency heating sections with heat utilization efficiency coefficients exceeding a preset threshold are selected from the heat utilization efficiency coefficients. Within the high-efficiency heating sections, the locations of heating sources are determined according to the spatial distribution of the antifreeze power demand. The spatial gradient change rate of the minimum antifreeze power demand is extracted from the spatial distribution of the antifreeze power demand. The locations where the spatial gradient change rate is greater than a preset change rate are identified as densely distributed heating source areas. The density of heating sources is increased within the densely distributed heating source areas. By combining the location of the heating sources in the high-efficiency heating section, the minimum antifreeze power requirement corresponding to the location of the heating sources, and the density of the heating sources in the densely distributed area, a dynamic spatial configuration scheme for antifreeze heave of heating sources in cold-region tunnels is generated, which includes the spatial location coordinates of the heating sources, the power configuration of the heating sources, and the spacing between the heating sources.

8. A dynamic heating system for preventing frost heave in cold-region tunnels based on multi-physics field coupling, used to implement the method as described in any one of claims 1-7, characterized in that, include: The first unit is used to obtain the temperature gradient, water content distribution, and freezing front advance speed at different depths in cold region tunnels. The second unit is used to perform coupled analysis of the temperature gradient and the moisture content distribution to obtain the moisture migration driving potential field. Based on the spatial distribution characteristics of the moisture migration driving potential field, the risk section of continuous migration and accumulation of moisture towards the freezing front in the cold region tunnel is identified, and a distribution map of the critical area of ​​frost heave risk is obtained. The third unit is used to establish a time priority sorting by calculating the dynamic equilibrium response time between the frozen front and the unfrozen water interface based on the freezing front advance speed of each segment in the distribution map of the critical area of ​​freezing heave risk, and to obtain a heating response sequence with time dimension constraints. The fourth unit is used to identify the dominant heat transfer path and heat dissipation path between the heating source and the freezing front based on the distribution law of the temperature gradient in the three-dimensional space of the cold region tunnel by tracking the evolution trajectory of the temperature isosurface, and to obtain the heat transfer path network topology. The fifth unit is used to calculate the heat utilization efficiency coefficient of different heating positions in each section by quantifying the contribution of temperature rise and energy attenuation ratio when a unit heating power is transferred to the freezing front along the current path in each path of the transmission path network topology. The sixth unit is used to generate a dynamic heating source spatial configuration scheme for cold-region tunnel frost heave prevention based on the heating response sequence and the heat utilization efficiency coefficient.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images