Deep circular cavern excavation response prediction method considering hot water force full coupling effect

Abstract

The invention provides a deep circular cavern excavation response prediction method considering a hot water force full coupling effect. The deep circular cavern excavation response prediction method comprises the steps that the state of an elastic-plastic solution interface is assumed; radial and tangential effective stresses at an elastic-plastic interface are obtained; performing layer-by-layer difference on the stress state at the elastic-plastic interface in the tunnel direction to obtain the stress of the ith layer of circular ring calculation unit, and obtaining the strain of the layer of circular ring calculation unit; updating to obtain calculation parameters such as damage variables of the (i + 1) th layer of circular ring calculation unit; if r is not greater than R0, continuing to repeat the above calculation process to the next layer of circular ring calculation unit; and if r is greater than or equal to R0, checking an error between a difference calculation value and a set value of the effective stress, the pore water pressure and the temperature at the excavation boundary, if the error is too large, updating a state assumed value at the elastic-plastic interface according to a gradient descent algorithm, and continuously performing inversion until the error accords with expectation to obtain a final solution. According to the method, the problem that an effective calculation model is lacked in the existing deep hot water force coupling environment cavern excavation design is solved.

Description

Technical Field

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

[0002] The use of underground space is gradually entering deep space, facing the thermal-hydraulic coupling environment. In the thermal-hydraulic coupling environment, the unloading effect caused by cavern excavation is different from that in the non-coupled environment. In general, the thermal-hydraulic coupling environment will aggravate the unloading effect of the surrounding rock, and the direct impact on the project is that the convergence deformation increases significantly and the surrounding rock fracture depth is greater. Originally, the prediction of convergence displacement and fracture zone radius after cavern excavation is particularly important for the construction of underground projects, and is the key parameter to guide cavern support. However, the evaluation of cavern excavation response under the thermal-hydraulic coupling environment still mainly relies on numerical simulation, and lacks an effective theoretical analysis model. It is undeniable that numerical simulation is a powerful tool for underground engineering construction. However, solving large-scale engineering problems, especially in the thermal-hydraulic coupling environment, often requires a long calculation time due to the nonlinearity of the problem. For such complex problems, the method of establishing a simplified model and solving it analytically is usually adopted. In fact, the preliminary design and most designs are based on relatively simple analysis models. Therefore, it is very necessary and urgent to establish an analysis model for underground cavern excavation that can effectively reflect the thermal-hydraulic coupling environment.

[0003] In the construction of deep underground projects, the problem of thermal-hydraulic coupling is often encountered. 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 geology, but the stress adjustment process of the surrounding rock around the underground opening is also significantly different from that under stress conditions alone. Accurately predicting the stability of the surrounding rock after cavern excavation is crucial to underground project construction and is a prerequisite for ensuring the safety of construction personnel and subsequent support design. Although many scholars have extended the analytical solution of tunnel mechanics, most of them are still in the context of non-coupled fields. However, as underground projects develop deeper, the surrounding rock is mostly in a coupled field environment. The experience of such engineering construction shows that the unloading effect of the surrounding rock will be intensified under the coupled environment. To this end, some scholars have approximated the solution by simply superimposing the various fields. However, the influence of water pressure and temperature on the unloading effect of the surrounding rock is not just a simple superposition of stress, permeability and thermal stress, but a real-time dynamic coupling between the three fields, that is, the interaction between stress, water pressure and temperature. Generally speaking, this dynamic coupling can be considered by establishing the dynamic relationship between rock damage and various physical quantities. Studies have shown that it is necessary to consider the coupling between surrounding rock damage and various fields in the prediction of surrounding rock behavior under the action of thermal-hydraulic coupling field. However, due to the interaction between various fields, it is very difficult to solve this series of problems. Therefore, the stability analysis of surrounding rocks under thermal-hydraulic coupling environment is currently mostly carried out using numerical simulation.

[0004] Obviously, numerical simulation is a powerful tool for solving underground engineering design and construction problems. However, numerical software usually takes a long time to solve large-scale engineering problems such as strongly nonlinear problems such as thermal-hydraulic coupling. Therefore, for such complex engineering problems, it is usually possible to establish a simplified model and solve it analytically, which can provide guidance for the subsequent preliminary support design. Summary of the invention

[0005] In view of 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 problem of the lack of an effective calculation model for the current deep thermal-hydraulic coupling environment cavern excavation design.

[0006] In order to achieve the above objectives, the technical solution adopted by the present invention is: a method for predicting the excavation response of a deep circular cavern considering the full coupling effect of thermal water and hydraulics, comprising the following steps:

[0007] S1, the state at the interface of the elastic-plastic solution is assumed 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. According to 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. According to the sum of the radial effective stress and the tangential effective stress and the yield criterion, the radial effective stress and the tangential effective stress at the elastic-plastic interface are calculated respectively;

[0011] S5, discretize the plastic region to obtain a number of equally spaced circular ring calculation units in the plastic region, a total of n+1 circular ring calculation nodes, and use the calculation result of S4 as the state value of the first circular ring calculation node;

[0012] S6. Calculate the state value of the (i+1)th layer of ring computing units according to the state value of the ring computing units in the i-th layer;

[0013] S7, judging whether the distance between the current calculation annular unit and the center of the cavern is greater than the excavation radius, if so, treating the i+1th layer annular calculation unit as a new i-th layer annular calculation unit, and returning to step S6, otherwise, entering 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 of 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 from the boundary to 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, and 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, 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 the tangential effective stress of the i-th layer of the annular calculation unit 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 the annular calculation unit, represents the radial effective stress of the i-th ring calculation unit, r i+1 represents the distance between the i+1th layer of 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 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, β i+1 represents the Biot coefficient of the i+1th layer of the 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 i-th ring calculation unit, represents the tangential effective stress of the i-layer annular calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the i-th ring computing unit, represents the friction angle of the i-th ring calculation unit, c represents the instantaneous cohesion of the surrounding rock;

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

[0025]

[0026] in, represents the tangential plastic strain of the i+1th layer of circular ring computing unit, represents the tangential plastic strain of the i-th ring computing unit, represents the radial plastic strain of the i-th ring calculation unit, represents the radial elastic strain of the i-th ring computing unit, 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 the annular computing unit, represents the tangential elastic strain of the i-th ring computing unit, represents the radial plastic strain of the i+1th layer of annular calculation unit, represents the radial plastic strain of the i-th layer of circular ring calculation unit, and ψ represents the shear dilatancy angle.

[0027] Further, the state value of the ring computing unit at the i+1th layer is calculated according to the state value of the ring computing unit at the ith 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] Combined with 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 units, the calculation parameters of the i+1th layer of circular ring calculation units are updated;

[0031] The original i+1th layer of circular ring calculation unit is regarded as the new i-th layer of circular ring calculation unit, and the state value of the new i+1th layer of circular ring calculation unit is continued until the distance between the circular ring calculation unit and the tunnel center is no greater than the tunnel excavation radius, and the calculation of the state value of the i+1th layer of circular ring calculation unit is completed, where the state value includes 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 ring computing unit and the tunnel center, k i represents the permeability of the i-th ring calculation unit;

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

[0036]

[0037] Among them, T i+1 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, Φ 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 the 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, 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] The present invention is intended to provide a theoretical analysis method for underground engineering construction in a thermal-hydraulic coupling environment. In view of the fact that the current analysis model for underground engineering construction in a thermal-hydraulic coupling environment fails to consider the full coupling effect between the three physical fields of heat, water and force, the present invention provides a calculation algorithm, which uses a gradient descent algorithm to invert the key parameters in the parameters, namely, water inflow, heat flux and plastic zone radius, thereby obtaining the calculation results considering the full coupling effect of heat, water and force. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The figure is a flow chart of the method of the present invention.

[0045] Figure 2 Schematic diagram of the elastoplastic 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 implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0049] Example

[0050] The present invention is intended to provide a theoretical analysis method for underground engineering construction in a thermal-hydraulic coupling environment. In view of the fact that the current analysis model for underground engineering construction in a thermal-hydraulic coupling environment fails to consider the full coupling effect between the three physical fields of heat, water and force, the present invention provides a calculation algorithm, which uses a gradient descent algorithm to invert the key parameters in the parameters, namely, water inflow, heat flux and plastic zone radius, thereby obtaining the calculation results considering the full coupling effect of heat, water and force.

[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 coupling effect of thermal water and hydraulics, and the implementation method thereof is as follows:

[0052] S1, the state at the interface of the elastic-plastic solution is assumed 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. According to 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. According to the sum of the radial effective stress and the tangential effective stress and the yield criterion, the radial effective stress and the tangential effective stress at the elastic-plastic interface are calculated respectively;

[0056] S5, discretize the plastic region to obtain a number of equally spaced circular ring calculation units in the plastic region, a total of n+1 circular ring calculation nodes, and use the calculation result of S4 as the state value of the first circular ring calculation node;

[0057] S6. Calculate the state value of the i+1th layer of ring computing units according to the state value of the i-th layer of ring computing units, 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] Combined with 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 units, the calculation parameters of the i+1th layer of circular ring calculation units are updated;

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

[0062] S7, judging whether the distance between the current calculation annular unit and the center of the cavern is greater than the excavation radius, if so, treating the i+1th layer annular calculation unit as a new i-th layer annular calculation unit, and returning to step S6, otherwise, entering 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 cave 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 radius of the plastic zone, R e It represents the calculation boundary of the model, that is, the radius of the elastic zone. It is considered that when the cavern is excavated under the thermal-hydraulic coupling 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 area and a plastic area. It is considered that the influence of the water pressure field, temperature field, and stress field stops at a sufficient 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 evolve dynamically only 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 the thermal-hydraulic coupling environment, in addition to the mechanical stress field, there are also thermal stress fields and seepage stress fields in the surrounding rock, which is the complexity of the thermal-hydraulic coupling environment problem. After the excavation of the cavern, the readjustment of the surrounding rock state in the non-coupled environment only needs to satisfy the balance of the mechanical stress field. Usually, the solution to this problem can be obtained by combining the Lame analytical solution in elastic mechanics and the yield criterion of the rock to obtain the stress state on the elastic-plastic interface. However, in the thermal-hydraulic coupling environment, due to the participation of thermal stress and seepage stress, the above solution ideas 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 are finally balanced independently of each other, and the three stress fields satisfy certain physical laws. This is the full thermal-hydraulic coupling effect. The states of the three physical fields can usually 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 decline of cohesion, based on the fact that the decline of rock strength after failure is mainly caused by the dissipation of cohesion, and the reduction of friction angle can be ignored. Therefore, 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, and 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, 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: elastic-plastic zone and 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, and 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 It represents the radius of the model calculation boundary, 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 radial effective stress, E represents elastic modulus, α represents thermal expansion coefficient, β represents 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 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 radial effective stress, and C5 and C6 represent stress distribution coefficient terms.

[0098] Substituting formula (11) into formula (6) we obtain:

[0099]

[0100] Think r = R e At this point, the influence of excavation disturbance is approximately 0. So substituting 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 geometry 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 equation (14) and equation (8) into equation (2), we can obtain:

[0106]

[0107] Solving for this yields:

[0108]

[0109] Among them, u represents the radial displacement, 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 circular 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 ring computing unit and the tunnel center, k i represents the permeability of the i-th ring computing unit, T i+1 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, Φ 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 i-th ring calculation unit, r i+1 represents the distance between the i+1th ring computing unit and the tunnel center, r irepresents the distance between the i-th ring computing unit and the tunnel center, represents the tangential effective stress of the calculation unit of the i-th ring, E represents the elastic modulus, α represents the thermal expansion coefficient, β represents the Biot coefficient, ΔT i+1 represents the temperature change of the i+1th ring computing unit, ΔT i Represents the temperature change of the i-th ring computing unit.

[0117] Combining the yield equation (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 the annular calculation unit, represents the radial effective stress of the i-th ring calculation unit, r i+1 represents the distance between the i+1th layer of 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 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, β i+1 represents the Biot coefficient of the i+1th layer of the 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 i-th ring calculation unit, represents the tangential effective stress of the i-layer annular calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the i-th ring computing unit, represents the friction angle of the i-th ring calculation unit, and c represents the instantaneous cohesion of the surrounding rock.

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

[0121]

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

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

[0124]

[0125] in, represents the tangential plastic strain of the i+1th layer of circular ring computing unit, represents the tangential plastic strain of the i-th ring computing unit, represents the radial plastic strain of the i-th ring calculation unit, represents the radial elastic strain of the i-th ring computing unit, 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 the annular computing unit, represents the tangential elastic strain of the i-th ring computing unit, represents the radial plastic strain of the i+1th layer of annular 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 are used to calculate the temperature, water pressure and stress of the i+1 layer ring by calculating the unit temperature, water pressure and stress of the i layer ring in the plastic zone. However, since the temperature, water pressure and stress at the elastic-plastic interface cannot be obtained directly, 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 radius of the plastic zone R p As unknown parameters in the model, all that needs to be done is to solve these three parameters using the boundary conditions in the model. 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 time, the temperature, stress, and water pressure on the elastic-plastic interface are obtained. By replacing the boundary conditions in Eq. 17 and Eq. 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 , we need to determine the correct initial conditions The calculated radial stress, temperature and pore pressure of the excavation boundary are matched with the actual boundary conditions. The correct initial conditions can be obtained by 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 of each iteration to ensure that the value of the objective function gradually decreases with each iteration. For this purpose, 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 represents the excavation boundary state parameter set according to the 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 hot-water-hydraulic coupling underground engineering. The 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 hot-water-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 the 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 calculation 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 project.

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, the state at the interface of the elastic-plastic solution is assumed 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. According to 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. According to the sum of the radial effective stress and the tangential effective stress and the yield criterion, the radial effective stress and the tangential effective stress at the elastic-plastic interface are calculated respectively; S5, discretize the plastic region to obtain a number of equally spaced circular ring calculation units in the plastic region, a total of n+1 circular ring calculation nodes, and use the calculation result of S4 as the state value of the first circular ring calculation node; S6. Calculate the state value of the (i+1)th layer of ring computing units according to the state value of the ring computing units in the i-th layer; S7, judging whether the distance between the current calculation annular unit and the center of the cavern is greater than the excavation radius, if so, treating the i+1th layer annular calculation unit as a new i-th layer annular calculation unit, and returning to step S6, otherwise, entering 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 of 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 from the boundary to 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 deep circular cavern excavation response prediction method considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The expression of the yield criterion is 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, and 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, 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 deep circular cavern excavation response prediction method 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 i-th layer annular calculation unit 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 the annular calculation unit, represents the radial effective stress of the i-th ring calculation unit, r i+1 represents the distance between the i+1th layer of 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 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, β i+1 represents the Biot coefficient of the i+1th layer of the 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 i-th ring calculation unit, represents the tangential effective stress of the i-layer annular calculation unit, represents the pore water pressure of the i-th ring calculation unit, T i represents the temperature of the i-th ring computing unit, represents the friction angle of the i-th ring calculation unit, c represents the instantaneous cohesion of the surrounding rock; The expressions of the tangential plastic strain and radial plastic strain of the i+1th layer of annular calculation unit are as follows: in, represents the tangential plastic strain of the i+1th layer of circular ring computing unit, represents the tangential plastic strain of the i-th ring computing unit, represents the radial plastic strain of the i-th ring calculation unit, represents the radial elastic strain of the i-th ring computing unit, 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 the annular computing unit, represents the tangential elastic strain of the i-th ring computing unit, represents the radial plastic strain of the i+1th layer of annular 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 deep circular cavern excavation response prediction method considering the full thermal-hydraulic coupling effect according to claim 1 is characterized in that: The state value of the ring computing unit at the i+1th layer is calculated according to the state value of the ring computing unit at the ith 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; Combined with 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 units, the calculation parameters of the i+1th layer of circular ring calculation units are updated; The original i+1th layer of circular ring calculation unit is regarded as the new i-th layer of circular ring calculation unit, and the state value of the new i+1th layer of circular ring calculation unit is continued until the distance between the circular ring calculation unit and the tunnel center is no greater than the tunnel excavation radius, and the calculation of the state value of the i+1th layer of circular ring calculation unit is completed, where the state value includes effective stress, plastic strain, pore water pressure and temperature.

6. The deep circular cavern excavation response prediction method 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 of 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 ring computing unit and the tunnel center, k i represents the permeability of the i-th ring calculation unit; The expression of the temperature of the i+1th layer ring calculation unit is as follows: Among them, T i+1 represents the temperature of the i+1th ring computing unit, T i represents the temperature of the i-th ring computing unit, Φ represents the heat flux, and λ i Represents the thermal conductivity of the i-th ring calculation unit.

7. 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 of the new elastic-plastic interface state assumed 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 the 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