Prediction method for deep circular cavern excavation response considering fully coupled thermal and hydraulic effects

Abstract

The present invention provides a method for predicting the excavation response of deep circular caverns that considers the full coupling of thermal and hydrothermal effects. The method includes assuming the state at the elastic-plastic solution interface; obtaining the radial and tangential effective stresses at the elastic-plastic interface; differentially determining the stress state at the elastic-plastic interface layer by layer toward the tunnel direction to obtain the stress of the i-th layer of circular ring calculation units, and obtaining the strain of the circular ring calculation units in this layer; updating the calculation parameters such as the damage variable of the i+1-th layer of circular ring calculation units; if r is not greater than R0, repeating the above calculation process for the next layer of circular ring calculation units; if r is greater than or equal to R0, checking the error between the differential calculated value and the set value of the effective stress, pore water pressure, and temperature at the excavation boundary. If the error is too large, updating the assumed value of the state at the elastic-plastic interface using a gradient descent algorithm, and continuing the inversion until the error meets the expected value to obtain the final solution. The present invention solves the current problem of the lack of an effective calculation model for the design of deep cavern excavation in thermal and hydrothermal coupled environments.

Description

Technical Field

[0001] The present invention belongs to the technical field of underground cavern construction, and in particular relates to a deep circular cavern excavation response prediction method considering full thermal-hydraulic coupling effects. Background Art

[0002] Underground space is increasingly being used at depth, facing a coupled hydrothermal and hydraulic environment. The unloading effect of cavern excavation in a coupled hydrothermal and hydraulic environment differs from that in an uncoupled environment. Generally speaking, a coupled hydrothermal and hydraulic environment exacerbates the unloading effect of the surrounding rock, directly impacting the project through a significant increase in convergence deformation and deeper fractures. Predicting convergence displacement and fracture zone radius after cavern excavation is crucial for underground engineering construction and serves as key parameters guiding cavern support. Currently, the assessment of cavern excavation response under coupled hydrothermal and hydraulic conditions still primarily relies on numerical simulation, lacking effective theoretical analytical models. While numerical simulation is undoubtedly a powerful tool for underground engineering construction, solving large-scale engineering problems, particularly in coupled hydrothermal and hydraulic environments, often requires significant computational time due to their nonlinear nature. For such complex problems, simplified models are typically developed and analytically solved. In fact, preliminary and most designs are based on relatively simple analytical models. Therefore, developing an analytical model that effectively reflects underground cavern excavation under coupled hydrothermal and hydraulic conditions is crucial and urgent for engineering.

[0003] Thermal-thermal coupling is a common problem encountered in deep underground engineering construction. Under the combined effects of stress, pore water pressure, and rock temperature, not only do the physical properties of the rock mass change relative to conventional geological conditions, but the stress adjustment process in the surrounding rock around underground openings also differs significantly from that under stress alone. Accurately predicting the stability of the surrounding rock after excavation is crucial for underground engineering construction, ensuring the safety of construction workers and the design of subsequent support structures. Although numerous researchers have extended analytical solutions to tunnel mechanics, most of these approaches are based on uncoupled fields. However, as underground engineering projects progress deeper, the surrounding rock often experiences coupled-field conditions. Experience with such projects indicates that the unloading effect of the surrounding rock is exacerbated in such a coupled environment. To this end, some researchers have employed approximate solutions through a simple superposition of various fields. However, the influence of water pressure and temperature on the unloading effect of the surrounding rock is not simply a superposition of stress, seepage force, and thermal stress, but rather a real-time dynamic coupling between these three fields: the interaction between stress, water pressure, and temperature. Generally speaking, this dynamic coupling can be considered by establishing a dynamic relationship between rock damage and various physical quantities. Research has shown that it is essential to consider the coupling between surrounding rock damage and various fields when predicting surrounding rock behavior under the action of thermal-hydraulic coupling fields. However, due to the interaction between the various fields, solving this problem is very difficult. Therefore, current analysis of surrounding rock stability under thermal-hydraulic coupling environments is mostly carried out using numerical simulation.

[0004] Numerical simulation is clearly a powerful tool for solving underground engineering design and construction problems. However, numerical software typically requires long runtimes to solve large-scale, highly nonlinear problems like thermal-hydraulic coupling. Therefore, for such complex engineering problems, simplified models can often be constructed and analytically solved to guide subsequent preliminary support design. Summary of the Invention

[0005] In response to the above-mentioned deficiencies in the prior art, the present invention provides a deep circular cavern excavation response prediction method considering the full thermal-hydraulic coupling effect, which solves the current problem of lack of an effective calculation model for deep thermal-hydraulic coupled environment cavern excavation design.

[0006] To achieve the above objectives, the present invention adopts a technical solution: a method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect, comprising the following steps:

[0007] S1. Assume the state at the interface of the elastic-plastic solution within the preset range of the deep circular cavern;

[0008] S2, calculating and obtaining the water pressure distribution function coefficient term C1 and the temperature distribution function coefficient term C3 in the elastic calculation area;

[0009] S3. Based on the assumed state at the elastic-plastic solution interface, the sum of the radial effective stress and the tangential effective stress at the elastic-plastic interface is calculated;

[0010] S4. Calculate the radial effective stress and tangential effective stress at the elastic-plastic interface based on the sum of the radial effective stress and the tangential effective stress and the yield criterion;

[0011] S5. Discretize the plastic region to obtain several equally spaced circular calculation units in the plastic region, with a total of n+1 circular calculation nodes. The calculation result of S4 is used as the state value of the first circular calculation node.

[0012] S6. Calculate the state value of the (i+1)th layer of ring computing units based on the state value of the i-th layer of ring computing units.

[0013] S7, determine whether the distance between the current calculation ring unit and the cavern center is greater than the excavation radius. If so, treat the (i+1)th layer ring calculation unit as a new (i)th layer ring calculation unit and return to step S6. Otherwise, proceed to step S8.

[0014] S8. Determine whether the relative errors between the pore water pressure, temperature, and radial effective stress obtained by differential calculation at the excavation boundary and the set conditions meet the requirements. If so, complete the response prediction of the deep circular cavern excavation. Otherwise, return to step S1 and use the gradient descent method to calculate and obtain a new assumed value of the elastic-plastic interface state.

[0015] Furthermore, the expression for the sum of the radial effective stress and the tangential effective stress is as follows:

[0016]

[0017] Among them, G(R p ) represents the sum of the radial effective stress and the tangential effective stress at the elastic-plastic interface, σ'0 represents the initial ground stress, R e Indicates the distance between the boundary and the center of the cavern, R p represents the distance between the elastic-plastic interface and the center of the cavern, E represents the elastic modulus of the surrounding rock, α represents the linear elastic density of the surrounding rock unit, C1 and C3 represent the temperature distribution coefficient term and the pore water pressure distribution coefficient term in the elastic region of the surrounding rock, v represents the Poisson's ratio, T p and T0 represent the ground temperature at the elastic-plastic interface and the initial ground temperature, respectively. w,p and P w,0 represent the pore water pressure at the elastic-plastic interface and the initial pore water pressure, respectively.

[0018] Furthermore, the yield criterion is expressed as follows:

[0019]

[0020] Where F represents the yield surface function, F(σ θ ,σ r )=0 means the surrounding rock unit reaches plasticity, σ θ and σ r denote the tangential effective stress and radial effective stress, respectively. represents the instantaneous friction angle, c represents the instantaneous cohesion of the surrounding rock, c ini represents the initial cohesion, γ p represents the equivalent plastic shear strain, n' represents the constant of the cohesion evolution curve, represents the tangential plastic strain, represents radial plastic strain, D represents damage variable, c res represents the residual cohesion, ω represents the instantaneous porosity, ω ini represents the initial porosity, ω res represents the residual porosity, k represents the instantaneous permeability, and k ini represents the initial permeability, λ represents the instantaneous thermal conductivity of the surrounding rock, and λ srepresents the thermal conductivity of the rock skeleton, λ w Represents the thermal conductivity of water.

[0021] Furthermore, the expressions of the radial effective stress and tangential effective stress of the annular calculation unit of the i-th layer are as follows:

[0022]

[0023] in, represents the radial effective stress of the i+1th layer of annular calculation unit, represents the tangential effective stress of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, r i+1 represents the distance between the i+1th layer circular calculation unit and the tunnel center, E represents the elastic modulus of the surrounding rock, β represents the Biot coefficient of the surrounding rock, α represents the thermal expansion coefficient of the surrounding rock, and r i represents the distance between the i-th ring computing unit and the tunnel center, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, β i+1 represents the Biot coefficient of the i+1th layer ring computing unit, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, represents the tangential effective stress of the i-layer ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the ring computing unit in the i-th layer, represents the friction angle of the i-th ring calculation unit, and c represents the instantaneous cohesion of the surrounding rock;

[0024] The expressions for the tangential plastic strain and radial plastic strain of the i+1th layer ring calculation unit are as follows:

[0025]

[0026] in, represents the tangential plastic strain of the i+1th layer of circular ring calculation unit, represents the tangential plastic strain of the i-th layer of ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, represents the radial elastic strain of the ring computing unit in the i-th layer, K ψ represents the plasticity multiplier, r i+1 represents the distance between the i+1th ring computing unit and the tunnel center, r iIndicates the distance between the i-th ring computing unit and the tunnel center, represents the tangential elastic strain of the i+1th layer of circular ring computing unit, represents the tangential elastic strain of the i-th layer of circular ring computing unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, and ψ represents the shear dilatancy angle.

[0027] Furthermore, the state value of the ring calculation unit in the i+1th layer is calculated based on the state value of the ring calculation unit in the i-th layer, which is specifically:

[0028] Based on the pore water pressure and temperature of the i-th layer of circular ring calculation unit, calculate the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit;

[0029] Combining the stress state and calculation parameters of the i-th layer of circular ring calculation unit and the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit, the effective stress and plastic strain of the i+1-th layer of circular ring calculation unit are obtained;

[0030] Using the plastic strain of the i+1th layer of circular ring calculation unit, the calculation parameters of the i+1th layer of circular ring calculation unit are updated;

[0031] The original i+1th layer of circular ring calculation units is regarded as the new i-th layer of circular ring calculation units, and the state values ​​of the new i+1th layer of circular ring calculation units are continued until the distance between the circular ring calculation unit and the tunnel center is no greater than the tunnel excavation radius. The calculation of the state values ​​of the i+1th layer of circular ring calculation units is completed, where the state values ​​include effective stress, plastic strain, pore water pressure, and temperature.

[0032] Furthermore, the expression of the pore water pressure of the i+1th layer of circular ring calculation unit is as follows:

[0033]

[0034] in, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, Q represents the seepage flux, r i represents the distance between the i-th ring computing unit and the tunnel center, r i+1 represents the distance between the i+1th layer ring computing unit and the tunnel center, k i represents the permeability of the ring calculation unit in the i-th layer;

[0035] The temperature of the i+1th ring calculation unit is expressed as follows:

[0036]

[0037] Among them, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, Φ represents the heat flux, and λ i Represents the thermal conductivity of the i-th ring calculation unit.

[0038] Furthermore, the expression of the new elastic-plastic interface state assumption value is as follows:

[0039]

[0040] LF(R p ,P w,p ,T p )=||(P s ,P w,s ,T s ) calu -(P s ,P w,s ,T s ) actu ||

[0041] Among them, (R p ,P w,p ,T) j represents the assumed solution of the last round of elastic-plastic interface state, (R p ,P w,p ,T) j+1 represents the new elastic-plastic interface state assumed solution, η represents the learning rate, Represents the gradient of the loss function, LF(R p ,P w,p ,T p ) represents the loss function, (P s ,P w,s ,T s ) calu Indicates the use of the assumed (R p ,P w,p ,T) The state of the excavation boundary obtained by calculation, (R p ,P w,p ,T) represents the assumed solution of elastic-plastic interface state, (P s ,P w,s ,T s ) actu Represents the excavation boundary state parameter set according to actual working conditions, R p represents the plastic zone radius, R e represents the radius of the elastic zone, T represents the temperature, T p represents the temperature of the elastic-plastic boundary, P s represents the radial support force at the excavation boundary, Pw,s represents the pore water pressure at the excavation boundary, T s Represents the temperature at the excavation boundary.

[0042] Beneficial effects of the present invention:

[0043] This invention aims to provide a theoretical analysis method for underground engineering construction in thermal-hydraulic coupled environments. Current analytical models for underground engineering construction in thermal-hydraulic coupled environments fail to account for the full coupling effects between the three physical fields of heat, water, and force. This invention proposes a calculation algorithm that uses a gradient descent algorithm to invert key parameters—namely, water inflow, heat flux, and the radius of the plastic zone—to obtain calculation results that account for the full thermal-hydraulic-force coupling effects. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 Flow chart of the method of the present invention.

[0045] Figure 2 Schematic diagram of the elastic-plastic analysis model for predicting the excavation response of a deep coupled environment cavern in this embodiment.

[0046] Figure 3 Schematic diagram of the fully coupled thermal water-hydraulic effect in this embodiment.

[0047] Figure 4 Schematic diagram of the discretization of the plastic zone in this embodiment. DETAILED DESCRIPTION

[0048] The specific embodiments of the present invention are described below to facilitate understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the appended claims, these changes are obvious, and all inventions and creations utilizing the concepts of the present invention are protected.

[0049] Example

[0050] This invention aims to provide a theoretical analysis method for underground engineering construction in thermal-hydraulic coupled environments. Current analytical models for underground engineering construction in thermal-hydraulic coupled environments fail to account for the full coupling effects between the three physical fields of heat, water, and force. This invention proposes a calculation algorithm that uses a gradient descent algorithm to invert key parameters—namely, water inflow, heat flux, and the radius of the plastic zone—to obtain calculation results that account for the full thermal-hydraulic-force coupling effects.

[0051] like Figure 1 As shown, the present invention provides a method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect, and its implementation method is as follows:

[0052] S1. Assume the state at the interface of the elastic-plastic solution within the preset range of the deep circular cavern;

[0053] S2, calculating and obtaining the water pressure distribution function coefficient term C1 and the temperature distribution function coefficient term C3 in the elastic calculation area;

[0054] S3. Based on the assumed state at the elastic-plastic solution interface, the sum of the radial effective stress and the tangential effective stress at the elastic-plastic interface is calculated;

[0055] S4. Calculate the radial effective stress and tangential effective stress at the elastic-plastic interface based on the sum of the radial effective stress and the tangential effective stress and the yield criterion;

[0056] S5. Discretize the plastic region to obtain several equally spaced circular calculation units in the plastic region, with a total of n+1 circular calculation nodes. The calculation result of S4 is used as the state value of the first circular calculation node.

[0057] S6. Calculate the state value of the ring computing unit in the (i+1)th layer according to the state value of the ring computing unit in the i-th layer, which is specifically:

[0058] Based on the pore water pressure and temperature of the i-th layer of circular ring calculation unit, calculate the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit;

[0059] Combining the stress state and calculation parameters of the i-th layer of circular ring calculation unit and the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit, the effective stress and plastic strain of the i+1-th layer of circular ring calculation unit are obtained;

[0060] Using the plastic strain of the i+1th layer of circular ring calculation unit, the calculation parameters of the i+1th layer of circular ring calculation unit are updated;

[0061] The original i+1th layer of circular ring calculation units is regarded as the new i-th layer of circular ring calculation units, and the state values ​​of the new i+1th layer of circular ring calculation units are continued until the distance between the circular ring calculation unit and the tunnel center is no greater than the tunnel excavation radius, completing the calculation of the state values ​​of the i+1th layer of circular ring calculation units, where the state values ​​include effective stress, plastic strain, pore water pressure, and temperature;

[0062] S7, determine whether the distance between the current calculation ring unit and the cavern center is greater than the excavation radius. If so, treat the (i+1)th layer ring calculation unit as a new (i)th layer ring calculation unit and return to step S6. Otherwise, proceed to step S8.

[0063] S8. Determine whether the relative errors between the pore water pressure, temperature, and radial effective stress obtained by differential calculation at the excavation boundary and the set conditions meet the requirements. If so, complete the response prediction of the deep circular cavern excavation. Otherwise, return to step S1 and use the gradient descent method to calculate and obtain a new assumed value of the elastic-plastic interface state.

[0064] In this embodiment, the present invention first establishes a cavern excavation response analysis model under a thermal-hydraulic coupling environment, such as Figure 2 As shown in the figure, T represents temperature, r represents the distance between the surrounding rock unit and the center of the cavern, and P w represents pore water pressure, σ'0 represents initial ground stress, P w,0 represents the initial pore water pressure, T0 represents the initial temperature, R0 represents the excavation radius of the cavern, and R p represents the plastic zone radius, R e Indicates the calculation boundary of the model, that is, the radius of the elastic zone. It is assumed that when the cavern is excavated under a thermal-hydraulic coupled environment, after the surrounding rock reaches a stable state, the area from the far end to the cavern can be divided into an elastic-plastic zone and a plastic zone. It is assumed that the influence of the water pressure field, temperature field, and stress field stops at a sufficiently far distance from the center of the cavern. This sufficiently far distance is marked as R in the figure. e The strength parameters (cohesion, friction angle) and coupling parameters (permeability, thermal conductivity) of the surrounding rock only evolve dynamically in the plastic zone, and the changes in the elastic-plastic zone are ignored.

[0065] Therefore, the above elastic-plastic analysis model under deep coupling environment is established to predict the surrounding rock response caused by cavern excavation. Figure 2 It can be seen that in a thermo-hydraulic coupled environment, the surrounding rock not only experiences mechanical stress fields, but also thermal and seepage stress fields, which underlies the complexity of the thermo-hydraulic coupled environment. After cavern excavation, readjustment of the surrounding rock state in an uncoupled environment only requires satisfying the equilibrium of the mechanical stress field. This problem is typically solved by combining the Lame analytical solution in elasticity with the yield criterion of rock to obtain the stress state at the elastic-plastic interface. However, in a thermo-hydraulic coupled environment, due to the involvement of thermal and seepage stresses, this approach cannot be achieved.

[0066] In fact, the mechanical stress field is closely related to the strength parameters of the rock (cohesion, friction angle), while the thermal stress field and seepage stress field are also related to the permeability and thermal conductivity of the surrounding rock, respectively. In the plastic zone, cohesion, friction angle, permeability, and thermal conductivity continue to evolve, so the mechanical stress field, thermal stress field, and seepage stress field evolve dynamically until they finally reach an independent balance among the mechanical stress field, thermal stress field, and seepage stress field, and the three stress fields satisfy certain physical laws. This is the full thermal-hydraulic coupling effect. Usually, the states of the three physical fields can be linked to each other through the damage variable D. The idea is as follows: Figure 3shown.

[0067] In this embodiment, the definition of the damage variable D varies from person to person. The present invention establishes an equation based on the degree of cohesion decay. This is based on the fact that the decay of rock strength after failure is primarily caused by the dissipation of cohesion, and the reduction in friction angle can be ignored. Thus, the overall result is:

[0068]

[0069] Where F represents the yield surface function, F(σ θ ,σ r )=0 means the surrounding rock unit reaches plasticity, σ θ and σ r denote the tangential effective stress and radial effective stress, respectively. represents the instantaneous friction angle, c represents the instantaneous cohesion of the surrounding rock, c ini represents the initial cohesion, γ p represents the equivalent plastic shear strain, n' represents the constant of the cohesion evolution curve, represents the tangential plastic strain, represents radial plastic strain, D represents damage variable, c res represents the residual cohesion, ω represents the instantaneous porosity, ω ini represents the initial porosity, ω res represents the residual porosity, k represents the instantaneous permeability, and k ini represents the initial permeability, λ represents the instantaneous thermal conductivity of the surrounding rock, and λ s represents the thermal conductivity of the rock skeleton, λ w Represents the thermal conductivity of water.

[0070] By Figure 2 The above equations are embedded in the computational model to consider the full coupling effect of thermal water and hydraulics. In addition, due to the dynamic evolution of mechanical parameters and coupling parameters of the surrounding rock in the plastic zone, the plastic zone of the computational model needs to be discretized, such as Figure 4 As shown. Therefore, we start to solve the model with the help of elastic-plastic theory. The specific solution process is divided into two parts: the elastic-plastic zone and the plastic zone:

[0071] 1) Elastic-plastic zone

[0072] The stress equilibrium equation of the surrounding rock is:

[0073]

[0074] Among them, σ' r represents the radial effective stress, σ' θrepresents the radial effective stress, r represents the distance between the surrounding rock unit and the center of the cavern, E represents the elastic modulus of the surrounding rock, ΔT represents the temperature variable, ΔT=T-T0, T represents the instantaneous temperature, T0 represents the initial temperature, β represents the Biot modulus, P w represents the pore water pressure, and α represents the thermal expansion coefficient of rock.

[0075] The distribution forms of temperature field and pore pressure field are:

[0076]

[0077] Among them, C1 and C2 represent the coefficient term and constant term of the temperature field distribution function in the elastic zone, respectively; C3 and C4 represent the coefficient term and constant term of the pore water pressure field distribution function, respectively.

[0078] According to the boundary conditions and It can be solved as:

[0079]

[0080] as well as:

[0081]

[0082] Among them, T p and T0 represent the temperature of the elastic-plastic boundary and the initial temperature respectively, P w,p and P w,0 denote the pore water pressure at the elastic-plastic boundary and the model calculation boundary, R p R represents the radius of the plastic zone, that is, the distance from the elastic-plastic boundary to the center of the cavern. e Indicates the model calculation boundary radius, that is, the distance from the model calculation boundary to the center of the cave.

[0083] Substituting equations (3) to (5) into equation (2), we obtain:

[0084]

[0085] Among them, σ' r represents the radial effective stress, σ' θ represents the radial effective stress, E represents the elastic modulus, α represents the thermal expansion coefficient, β represents the Biot coefficient, and r represents the distance between the surrounding rock unit and the center of the cavern.

[0086] On the other hand, the strain compatibility equation in the plane strain problem is:

[0087]

[0088] Among them, ε θ represents the tangential plastic strain, ε rrepresents the radial plastic strain, and r represents the distance between the surrounding rock unit and the center of the cavern.

[0089] The relationship between strain and stress satisfies Hooke's law:

[0090]

[0091] Where v is Poisson's ratio. Substituting Equation 8 into Equation 7, we can obtain:

[0092]

[0093] Substituting formula (9) into formula (2), we can obtain:

[0094]

[0095] Solving the above equations, we can obtain:

[0096]

[0097] Among them, σ' r represents the radial effective stress, σ' θ represents the radial effective stress, and C5 and C6 represent the stress distribution coefficient terms.

[0098] Substituting formula (11) into formula (6) yields:

[0099]

[0100] Think r = R e At this point, the influence of excavation disturbance is approximately 0. So we substitute the boundary conditions of the elastic-plastic zone into then:

[0101]

[0102] In order to determine the displacement of the surrounding rock in the elastic-plastic zone, it is necessary to use the unit geometry equation, Hooke's equation and stress equilibrium equation. The geometric equation of the rock unit is:

[0103]

[0104] Among them, ε r represents radial strain, ε θ represents the tangential strain, u represents the radial displacement, and r represents the distance between the surrounding rock unit and the center of the cavern.

[0105] Substituting equations (14) and (8) into equation (2), we can obtain:

[0106]

[0107] The solution is:

[0108]

[0109] Where u represents the radial displacement, and C7 and C8 represent the coefficient terms of the elastic zone displacement distribution function.

[0110] 2) Plastic zone

[0111] After entering the plastic zone, the mechanical parameters and coupling parameters of the rock mass at different distances from the tunnel are constantly evolving, so it is impossible to give a closed solution. Therefore, the plastic zone is discretized into n ring units. When n is large enough, it is assumed that the mechanical parameters and coupling parameters in each ring do not change, such as Figure 4 As shown. Therefore, within the plastic zone, the pore water pressure and rock temperature of the i+1th ring calculation unit can be obtained by integrating the pore water pressure and rock temperature of the i-th ring calculation unit through Darcy's law and Fourier's law respectively:

[0112]

[0113] in, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, Q represents the seepage flux, r i represents the distance between the i-th ring computing unit and the tunnel center, r i+1 represents the distance between the i+1th layer ring computing unit and the tunnel center, k i represents the permeability of the i-th ring calculation unit, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, Φ represents the heat flux, and λ i Represents the thermal conductivity of the i-th ring calculation unit.

[0114] The calculation of effective stress is relatively complicated. In the plastic zone, that is, R0 <r<R p In the interval, since n ring units with very small thickness are divided, it can be approximately considered that the coupling parameters do not change in the ring. At this time, the stress balance equation (2) of the surrounding rock unit can be rewritten in differential form:

[0115]

[0116] in, represents the radial effective stress of the i+1th layer of annular calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, r i+1 represents the distance between the i+1th ring computing unit and the tunnel center, r iIndicates the distance between the i-th ring computing unit and the tunnel center, represents the effective stress of the tangential direction of the ring calculation unit in the i-th layer, E represents the elastic modulus, α represents the thermal expansion coefficient, β represents the Biot coefficient, ΔT i+1 Indicates the temperature change of the i+1th ring computing unit, ΔT i Indicates the temperature change of the ring computing unit in the i-th layer.

[0117] Combining the yield equations (1), we get:

[0118]

[0119] in, represents the radial effective stress of the i+1th layer of annular calculation unit, represents the tangential effective stress of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, r i+1 represents the distance between the i+1th layer circular ring calculation unit and the tunnel center, E represents the elastic modulus of the surrounding rock, β represents the Biot coefficient of the surrounding rock, α represents the thermal expansion coefficient of the surrounding rock, and r i represents the distance between the i-th ring computing unit and the tunnel center, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, β i+1 represents the Biot coefficient of the i+1th layer ring computing unit, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, represents the tangential effective stress of the i-layer ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the ring computing unit in the i-th layer, represents the friction angle of the i-th layer circular ring calculation unit, and c represents the instantaneous cohesion of the surrounding rock.

[0120] Similarly, the differential form of the compatibility equation (7) is:

[0121]

[0122] in, represents the tangential plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the tangential plastic strain of the i-th layer of ring calculation unit, represents the tangential elastic strain of the i-th layer of circular ring computing unit, represents the radial elastic strain of the ring computing unit in the i-th layer, K ψ represents the plasticity multiplier.

[0123] Therefore, when the stress of the ring calculation unit in the i-th layer is known, the plastic strain calculation equation of the ring calculation unit in the i+1 layer can be further obtained:

[0124]

[0125] in, represents the tangential plastic strain of the i+1th layer of circular ring calculation unit, represents the tangential plastic strain of the i-th layer of ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, represents the radial elastic strain of the ring computing unit in the i-th layer, K ψ represents the plasticity multiplier, r i+1 represents the distance between the i+1th ring computing unit and the tunnel center, r i Indicates the distance between the i-th ring computing unit and the tunnel center, represents the tangential elastic strain of the i+1th layer of circular ring computing unit, represents the tangential elastic strain of the i-th layer of circular ring computing unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, and ψ represents the shear dilatancy angle.

[0126] The above formulas provide the calculation formula for calculating the temperature, water pressure, and stress of the i+1 ring layer by calculating the temperature, water pressure, and stress of the i-th ring layer in the plastic zone. However, because the temperature, water pressure, and stress at the elastic-plastic interface cannot be directly obtained, some algorithms are needed to obtain the final solution of the model.

[0127] In order to obtain the solution of the model, the temperature T at the elastic-plastic interface is p , pore water pressure P w,p And the plastic zone radius R p As unknown parameters in the model, all that needs to be done is to use the boundary conditions in the model to solve these three parameters. By combining equations (11) and (12), it can be seen that the sum of the radial stress and the tangential stress in the elastic-plastic zone is a function of the distance from this point to the center of the chamber, that is:

[0128]

[0129] Where C6 represents the coefficient term of the stress distribution function, which is calculated using formula (13).

[0130] So we can get r = R p Department,

[0131] By combining the yield equation (1), we can obtain the effective stress at the elastic-plastic interface: At this point, the temperature, stress, and water pressure at the elastic-plastic interface are obtained. By replacing the boundary conditions in Equations 17 and 18 with the boundary conditions of the model's elastic-plastic region, the seepage rate Q and heat flux Φ in the model can be calculated.

[0132] It can be seen that through the above calculation process, a set of initial conditions Corresponding to a set of excavation boundary states (T, σ r ,P w ) r＝R0 , the correct initial conditions need to be determined The calculated radial stress, temperature, and pore pressure at the excavation boundary are matched to the actual boundary conditions. The correct initial conditions can be obtained using the gradient descent method. The gradient descent method is an iterative optimization algorithm used to find the minimum value of a function. It uses the direction of the negative gradient of the loss function to determine the search direction for each iteration, ensuring that the value of the objective function gradually decreases with each iteration. To this end, the loss function can be defined as:

[0133] LF(R p ,P w,p ,T p )=||(P s ,P w,s ,T s ) calu -(P s ,P w,s ,T s ) actu || (23)

[0134] Therefore, after each iterative calculation, the calculation parameters at the new elastic-plastic interface are:

[0135]

[0136] Among them, (R p ,P w,p ,T) j represents the assumed solution of the last round of elastic-plastic interface state, (R p ,P w,p ,T) j+1 represents the new elastic-plastic interface state assumed solution, η represents the learning rate, Represents the gradient of the loss function, LF(R p ,P w,p ,T p ) represents the loss function, (P s ,P w,s ,T s ) calu Indicates the use of the assumed (R p ,P w,p ,T) The state of the excavation boundary obtained by calculation, (R p ,P w,p ,T) represents the assumed solution of elastic-plastic interface state, (P s ,P w,s ,T s ) actu Indicates the excavation boundary state parameters set according to actual working conditions, R p represents the radius of the plastic zone, R e represents the radius of the elastic zone, T represents the temperature, T p represents the temperature of the elastic-plastic boundary, P s represents the radial effective stress at the excavation boundary, i.e. the radial support stress, P w,s represents the pore water pressure at the excavation boundary, T s Represents the temperature at the excavation boundary.

[0137] Finally, the solution of the elastic-plastic solution model considering the full coupling effect of thermal water and hydraulics was obtained.

[0138] In summary, the present invention provides a theoretical analysis method for the construction of hydrothermal-hydraulic coupling underground engineering. Key engineering parameters such as the plastic zone radius and the convergence displacement of the cavern wall can be obtained through model calculation, providing parameter support for the support design of underground engineering construction. More importantly, studies have shown that if the full coupling effect is not considered in hydrothermal-hydraulic coupling engineering, it will be dangerous for the engineering. The reason is that if the full coupling effect is ignored, the unloading effect of the surrounding rock will be weakened, resulting in a reduction in key support parameters. For example, the plastic zone radius is often used as a key design parameter for the anchor length in anchor support, and ignoring the full coupling effect will cause the calculated result of the plastic zone radius to be too small, resulting in the anchor installed during the cavern excavation process to be too short, thereby causing support failure, delaying the construction period, and even causing casualties. The calculation model proposed in the present invention takes into account the full coupling effect, so the calculation result is relatively safe for the engineering.

Claims

1. A method for predicting the excavation response of a deep circular cavern considering the full coupling effect of thermal water and hydraulics, characterized in that: The following steps are involved: S1. Assume the state at the interface of the elastic-plastic solution within the preset range of the deep circular cavern; S2, calculating and obtaining the water pressure distribution function coefficient term C1 and the temperature distribution function coefficient term C3 in the elastic calculation area; S3. Based on the assumed state at the elastic-plastic solution interface, the sum of the radial effective stress and the tangential effective stress at the elastic-plastic interface is calculated; S4. Calculate the radial effective stress and tangential effective stress at the elastic-plastic interface based on the sum of the radial effective stress and the tangential effective stress and the yield criterion; S5. Discretize the plastic region to obtain several equally spaced circular calculation units in the plastic region, with a total of n+1 circular calculation nodes. The calculation result of S4 is used as the state value of the first circular calculation node. S6. Calculate the state value of the (i+1)th layer of ring computing units based on the state value of the i-th layer of ring computing units. S7, determine whether the distance between the current calculation ring unit and the cavern center is greater than the excavation radius. If so, treat the (i+1)th layer ring calculation unit as a new (i)th layer ring calculation unit and return to step S6. Otherwise, proceed to step S8. S8. Determine whether the relative errors between the pore water pressure, temperature, and radial effective stress obtained by differential calculation at the excavation boundary and the set conditions meet the requirements. If so, complete the response prediction of the deep circular cavern excavation. Otherwise, return to step S1 and use the gradient descent method to calculate and obtain a new assumed value of the elastic-plastic interface state.

2. The deep circular cavern excavation response prediction method considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The expression for the sum of the radial effective stress and the tangential effective stress is as follows: Among them, G(R p ) represents the sum of the radial effective stress and the tangential effective stress at the elastic-plastic interface, σ'0 represents the initial ground stress, R e Indicates the distance between the boundary and the center of the cavern, R p represents the distance between the elastic-plastic interface and the center of the cavern, E represents the elastic modulus of the surrounding rock, α represents the linear elastic density of the surrounding rock unit, C1 and C3 represent the temperature distribution coefficient term and the pore water pressure distribution coefficient term in the elastic region of the surrounding rock, v represents the Poisson's ratio, T p and T0 represent the ground temperature at the elastic-plastic interface and the initial ground temperature, respectively. w,p and P w,0 represent the pore water pressure at the elastic-plastic interface and the initial pore water pressure, respectively.

3. The method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The yield criterion is expressed as follows: Where F represents the yield surface function, F(σ θ ,σ r )=0 means the surrounding rock unit reaches plasticity, σ θ and σ r denote the tangential effective stress and radial effective stress, respectively. represents the instantaneous friction angle, c represents the instantaneous cohesion of the surrounding rock, c ini represents the initial cohesion, γ p represents the equivalent plastic shear strain, n' represents the constant of the cohesion evolution curve, represents the tangential plastic strain, represents radial plastic strain, D represents damage variable, c res represents the residual cohesion, ω represents the instantaneous porosity, ω ini represents the initial porosity, ω res represents the residual porosity, k represents the instantaneous permeability, and k ini represents the initial permeability, λ represents the instantaneous thermal conductivity of the surrounding rock, and λ s represents the thermal conductivity of the rock skeleton, λ w Represents the thermal conductivity of water.

4. The method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The expressions of the radial effective stress and tangential effective stress of the ring calculation unit in the i-th layer are as follows: in, represents the radial effective stress of the i+1th layer of annular calculation unit, represents the tangential effective stress of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, r i+1 represents the distance between the i+1th layer circular ring calculation unit and the tunnel center, E represents the elastic modulus of the surrounding rock, β represents the Biot coefficient of the surrounding rock, α represents the thermal expansion coefficient of the surrounding rock, and r i represents the distance between the i-th ring computing unit and the tunnel center, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, β i+1 represents the Biot coefficient of the i+1th layer ring computing unit, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the radial effective stress of the ring calculation unit in the i-th layer, represents the tangential effective stress of the i-layer ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the ring computing unit in the i-th layer, represents the friction angle of the i-th layer of circular ring calculation unit, and c represents the instantaneous cohesion of the surrounding rock; The expressions for the tangential plastic strain and radial plastic strain of the i+1th layer ring calculation unit are as follows: in, represents the tangential plastic strain of the i+1th layer of circular ring calculation unit, represents the tangential plastic strain of the i-th layer of ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, represents the radial elastic strain of the ring computing unit in the i-th layer, K ψ represents the plasticity multiplier, r i+1 represents the distance between the i+1th ring computing unit and the tunnel center, r i Indicates the distance between the i-th ring computing unit and the tunnel center, represents the tangential elastic strain of the i+1th layer of circular ring computing unit, represents the tangential elastic strain of the i-th layer of circular ring computing unit, represents the radial plastic strain of the i+1th layer of circular ring calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, and ψ represents the shear dilatancy angle.

5. The method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The state value of the ring computing unit in the i+1th layer is calculated based on the state value of the ring computing unit in the i-th layer, which is specifically: Based on the pore water pressure and temperature of the i-th layer of circular ring calculation unit, calculate the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit; Combining the stress state and calculation parameters of the i-th layer of circular ring calculation unit and the pore water pressure and temperature of the i+1-th layer of circular ring calculation unit, the effective stress and plastic strain of the i+1-th layer of circular ring calculation unit are obtained; Using the plastic strain of the i+1th layer of circular ring calculation unit, the calculation parameters of the i+1th layer of circular ring calculation unit are updated; The original i+1th layer of circular ring calculation units is regarded as the new i-th layer of circular ring calculation units, and the state values ​​of the new i+1th layer of circular ring calculation units are continued until the distance between the circular ring calculation unit and the tunnel center is no greater than the tunnel excavation radius. The calculation of the state values ​​of the i+1th layer of circular ring calculation units is completed, where the state values ​​include effective stress, plastic strain, pore water pressure, and temperature.

6. The method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect according to claim 5 is characterized in that: The expression of the pore water pressure of the i+1th layer circular ring calculation unit is as follows: in, represents the pore water pressure of the i+1th layer of circular ring calculation unit, represents the pore water pressure of the i-th ring calculation unit, Q represents the seepage flux, r i represents the distance between the i-th ring computing unit and the tunnel center, r i+1 represents the distance between the i+1th layer ring computing unit and the tunnel center, k i represents the permeability of the ring calculation unit in the i-th layer; The temperature of the i+1th ring calculation unit is expressed as follows: Among them, T i+1 Indicates the temperature of the i+1th ring computing unit, T i represents the temperature of the ring computing unit in the i-th layer, Φ represents the heat flux, and λ i Represents the thermal conductivity of the i-th ring calculation unit.

7. The method for predicting the excavation response of a deep circular cavern considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The expression of the new elastic-plastic interface state assumption value is as follows: LF(R p ,P w,p ,T p )=||(P s ,P w,s ,T s ) calu -(P s ,P w,s ,T s ) actu || Among them, (R p ,P w,p ,T) j represents the assumed solution of the last round of elastic-plastic interface state, (R p ,P w,p ,T) j+1 represents the new elastic-plastic interface state assumed solution, η represents the learning rate, Represents the gradient of the loss function, LF(R p ,P w,p ,T p ) represents the loss function, (P s ,P w,s ,T s ) calu Indicates the use of the assumed (R p ,P w,p ,T) The state of the excavation boundary obtained by calculation, (R p ,P w,p ,T) represents the assumed solution of elastic-plastic interface state, (P s ,P w,s ,T s ) actu Represents the excavation boundary state parameter set according to actual working conditions, R p represents the radius of the plastic zone, R e represents the radius of the elastic zone, T represents the temperature, T p represents the temperature of the elastic-plastic boundary, P s represents the radial support force at the excavation boundary, P w,s represents the pore water pressure at the excavation boundary, T s Represents the temperature at the excavation boundary.

Images