A numerical method for stress-seepage coupling in single joints

Abstract

The present invention discloses a single-joint stress-seepage coupling numerical calculation method, which belongs to the field of fluid-solid coupling numerical calculation, including: obtaining the input parameters of the single-joint stress-seepage coupling, discretizing the input parameters, and obtaining discretized units, wherein the discretized units include: contact units and non-contact units; based on the contact units, establishing contact equations, and based on the contact equations, obtaining contact stress and deformation; based on the deformation, updating the opening of the non-contact units to obtain updated results, performing contact judgment on the updated results, establishing the Reynolds equation of the non-contact units based on the judgment results, and obtaining the pressure and flow rate of the non-contact units based on the Reynolds equations, and judging the convergence based on the pressure. If it does not converge, continue to establish the contact equation until it converges. The present invention can solve the technical problems of high mesh requirements and long solution time, has fast calculation speed, simple mesh division, and realizes computer automatic processing.

Description

Technical Field

[0001] The invention belongs to the field of fluid-solid coupling numerical calculation, and in particular relates to a single-joint stress-seepage coupling numerical calculation method. Background Art

[0002] There is a coupling effect between the seepage field and the stress field of the rock mass. While the displacement field and the stress field are affected by the seepage load, they also affect the distribution of the seepage head in the rock mass. This is because there are many joints distributed in the rock mass. The stress action causes changes in the geometric characteristics of the joints, such as the opening or closing of the joint opening, the change of the joint roughness, etc., which affects the permeability of the joints and causes a significant change in the permeability of the entire rock mass. Changes in the seepage field will change the seepage force, thereby affecting the joint stress field, and the two affect each other. Joint seepage and stress coupling are important and hot topics in rock hydraulics.

[0003] Due to the roughness of the rock joint surface, the stress-seepage coupling mechanism of the joint is very complicated. Many experts and scholars at home and abroad generalize the joint roughness into one or two variables to establish a stress-seepage coupling model. These models ignore the influence of the microscopic protrusion distribution and the opening distribution of the joint surface on the stress-seepage coupling mechanism, which leads to inevitable errors. In order to accurately capture the influence of joint roughness, many scholars use numerical techniques such as CFD and discrete element method to establish a refined model of the joint and simulate the joint stress-seepage coupling process. However, this numerical simulation method requires a lot of time to establish a discrete grid model of the joint, and the numerical solution has high requirements on the grid, the solution time is long, and it often encounters non-convergence problems. Therefore, CFD simulation is more troublesome to apply. Summary of the invention

[0004] The purpose of the present invention is to provide a single joint stress-seepage coupling numerical calculation method to solve the problems of high mesh requirements and long solution time in the above-mentioned prior art.

[0005] To achieve the above object, the present invention provides a single joint stress-seepage coupling numerical calculation method, comprising the following steps:

[0006] Obtaining input parameters of single-joint stress-seepage coupling, and discretizing the input parameters to obtain discretized units, wherein the discretized units include: contact units and non-contact units;

[0007] Based on the contact unit, a contact equation is established, and based on the contact equation, contact stress and deformation are obtained;

[0008] Based on the deformation, the opening of the non-contact unit is updated to obtain an updated result, and a contact judgment is performed on the updated result. The Reynolds equation of the non-contact unit is established based on the judgment result. Based on the Reynolds equation, the pressure and flow rate of the non-contact unit are obtained. The convergence is judged based on the pressure. If it does not converge, the contact equation continues to be established until it converges.

[0009] Preferably, the process of discretizing the input parameters includes:

[0010] The input parameters of the single joint stress-seepage coupling are obtained, wherein the input parameters include: upper joint coordinates, lower joint coordinates, elevation, initial opening and contact radius, and the upper joint and the lower joint are discretized by using rectangular units to obtain discretized units.

[0011] Preferably, the contact equation includes: a force balance equation, a deformation coordination equation and an elasticity equation.

[0012] Preferably, before establishing the contact equation, the following steps are also included:

[0013] An initialization judgment is performed on the discretization units to determine whether the discretization units are in contact. If so, they are contact units; if not, they are non-contact units.

[0014] Preferably, the process of obtaining the contact stress and deformation includes:

[0015] Based on the force balance equation, the Gaussian iteration method is used to solve the force balance equation to obtain the contact stress and deformation of the protrusion.

[0016] Preferably, the process of performing contact judgment on the update result includes:

[0017] If the opening degree of the update result is less than 0, it is a contact unit; if the opening degree of the update result is greater than 0, it is a non-contact unit.

[0018] Preferably, the process of obtaining the pressure and flow rate of the uncontacted unit comprises:

[0019] Based on the Reynolds equation and the cubic law, the Galerkin method is used to establish an equivalent weak integral form of the Reynolds equation, and a shape function is established using a rectangular unit with bidirectional node interpolation. Based on the shape function and the equivalent weak integral form, the discrete form of each rectangular unit and a set of equations for all unknown node water pressures are obtained. Based on the set of equations, the pressure and flow rate of the uncontacted units are obtained.

[0020] Preferably, the process of judging convergence based on the pressure includes:

[0021] If the opening and pressure of the non-contacting unit satisfy the convergence relationship, it is converged; if the opening and pressure of the non-contacting unit do not satisfy the convergence relationship, it is not converged.

[0022] The technical effects of the present invention are:

[0023] The present invention provides a single-joint stress-seepage coupling numerical calculation method, which obtains input parameters of the single-joint stress-seepage coupling, discretizes the input parameters to obtain discretized units, wherein the discretized units include: contact units and non-contact units; based on the contact units, a contact equation is established, and based on the contact equation, contact stress and deformation are obtained; based on the deformation, the opening of the non-contact unit is updated to obtain an updated result, contact judgment is performed on the updated result, and a Reynolds equation of the non-contact unit is established based on the judgment result, and based on the Reynolds equation, the pressure and flow rate of the non-contact unit are obtained, and convergence is judged based on the pressure. If it does not converge, the contact equation continues to be established until convergence.

[0024] The present invention can solve the technical problems of high mesh requirements, long solution time and non-convergence in the prior art. The present invention performs stress-seepage coupling simulation on a certain crack to obtain an opening distribution diagram, a flow velocity diagram and a pressure diagram under different normal stresses. The present invention has a fast calculation speed, simple mesh division, fast division speed, and realizes computer automatic processing. The present invention takes into account the influence of crack contact on seepage, making the calculation more reasonable. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0026] Figure 1 is a flow chart of a calculation method in an embodiment of the present invention;

[0027] Figure 2 Schematic diagram of the contact area and non-contact area of ​​the discrete crack unit in the embodiment of the present invention;

[0028] Figure 3 A schematic diagram of cracks obtained by splitting a rock block in an embodiment of the present invention;

[0029] Figure 4 Schematic diagram of fracture morphology, elevation and opening measurement in an embodiment of the present invention;

[0030] Figure 5 Schematic diagram of the opening distribution under no stress in an embodiment of the present invention;

[0031] Figure 6Flow velocity diagram under different stresses in an embodiment of the present invention;

[0032] Figure 7 Graph showing pressure distribution under different stresses in an embodiment of the present invention. DETAILED DESCRIPTION

[0033] It should be noted that, in the absence of conflict, the embodiments and features in the embodiments of the present application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0034] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0035] Embodiment 1

[0036] like Figure 1 As shown, this embodiment provides a single joint stress-seepage coupling numerical calculation method, including:

[0037] S1. Input parameters: input upper and lower joint coordinates (x, y) and height h (x, y) respectively

[0038] S2. Discretization: Use rectangular units to discretize the upper and lower joints, and each unit is recorded as (i, j), such as Figure 2 shown.

[0039] S3. Initialization judgment: judge whether each unit is in contact, and the contact unit is recorded as (i, j) m , the uncontacted unit is recorded as (i,j) n

[0040] S4. Establish contact area (i, j) m The force balance equation, deformation coordination equation and elastic equation are solved to obtain the magnitude of contact unit force and displacement.

[0041] The force balance equation is as follows:

[0042]

[0043] Where F(i,j) m is the stress distribution, P(i,j) n is the pore water pressure distribution, R1(i,j) is the upper joint radius, R2(i,j) is the lower joint radius, see Figure 2 , R is the equivalent radius of the protrusion, R = (R1(i,j)+R2(i,j)) / 2, Δu is the deformation, and v is the flow velocity.

[0044] For the force balance equations, the Gaussian iteration method is used to solve the stress distribution F(i,j) of the protrusion. m and deformation Δu.

[0045] S5. According to the deformation value Δu obtained in the previous step, use the opening calculation formula to calculate the opening b(i,j) of the non-contact area n to update.

[0046] The calculation formula of opening is as follows:

[0047] b(i,j) n =b(i,j) n -Δu k (2)

[0048] S6. Determine whether each unit is in contact based on the calculated results: if b(i,j) n <0, the unit is in contact, recorded as (i,j) m , the opening is 0; if b(i,j) n > 0, the unit is not in contact, and the uncontacted unit is recorded as (i, j) n

[0049] S7. Create uncontacted element (i, j) n The overall conduction equation is discretized from the Reynolds equation and solved to obtain the pressure and flow velocity in the uncontacted area.

[0050] Considering that the water flow satisfies the Reynolds equation and the cubic law, the Reynolds equation and the cubic law are:

[0051]

[0052]

[0053] The Galerkin method is used to establish the equivalent weak integral form of the Reynolds equation:

[0054]

[0055] The shape functions are created using rectangular elements with bidirectional node interpolation:

[0056]

[0057] Substituting equation (6) into equation (5), we get the discrete form of each rectangular unit:

[0058] K e P e =F e (7)

[0059]

[0060] By adding the coefficients of the corresponding numbers in the unit matrix, the total conduction matrix in equation (5) can be formed, and the equation system for solving the water pressure of all unknown nodes can be obtained:

[0061] KP=F (11)

[0062] Among them, the elements of [K] Constant term The preconditioned conjugate gradient method is used to solve equation (11) to obtain the water pressure on the uncontacted area unit:

[0063] S8. Solve for the pressure value P(i,j) of the uncontacted unit m

[0064] S9. Determine whether it converges: If || b k -b k+1 ||<1e-4, ||P k -P k+1 ||<1e-4 converges and the calculation is completed, where b k and b k+1 They represent the opening size of the k and k+1 iteration steps of the non-contact unit, P k and P k+1 They represent the pressure magnitudes of the k and k+1 iteration steps of the non-contact unit respectively;

[0065] If it does not converge, update F(i,j) of the contacting and non-contacting elements. m and P(i,j) n , establish the force balance equation, deformation coordination equation and elastic equation of the contact unit again and solve them to obtain the contact unit force and deformation, and repeat S5-S8 until convergence.

[0066] S10. End calculation.

[0067] The specific implementation cases are as follows:

[0068] ①First, the 150mm×150mm×150mm granite specimen was Brazilian split to obtain fresh and complete cracks, such as Figure 3 As shown;

[0069] ② Crack morphology data measurement, initial opening and contact radius measurement. First, stick marking points (such as Figure 4 As shown in the figure, the rough surface of the upper and lower cracks and the contour under the "closed" state are scanned by a 3D optical scanning system to obtain point cloud data; then the data is imported into the Geomaic studio software, and the "stitching" function of the software is used to match the upper and lower crack marking points to "simulate" the closed state of the crack ( Figure 4); finally, take out the upper and lower crack surfaces, calculate the upper and lower morphological elevations, and then subtract the upper and lower joint elevations to obtain the initial opening distribution data of the closed crack. If the upper joint minus the lower joint at a certain point is less than 0, it means contact, and the opening at this point is 0. The lowest elevation of the joint above or below the contact radius is obtained by subtracting the lowest elevation at each point to obtain the contact radius value.

[0070] ③ Input parameters and discretization: Input the upper and lower joint coordinates (x, y), elevation h(x, y), initial opening b(i, j), contact radius R1(x, y) and R2(x, y) respectively. Use rectangular units to discretize the upper and lower joints, and each unit is recorded as (i, j), such as Figure 2 shown.

[0071] ④ Initialization judgment: judge whether each unit is in contact, and the contact unit is recorded as (i, j) m , the uncontacted unit is recorded as (i,j) n

[0072] ⑤Establish contact area (i,j) m The force balance equation, deformation coordination equation and elastic equation are solved to obtain the magnitude of contact unit force and displacement. At the same time, the new opening value is updated and whether there is contact is determined.

[0073] ⑥Establish uncontacted unit (i,j) n The overall conduction equation of the Reynolds equation is discretized and solved to obtain the pressure and flow velocity of the uncontacted area. The pressure value of the uncontacted unit is solved.

[0074] ⑦ Determine whether it converges: If ||b k -b k+1 ||<1e-4, ||P k -P k+1 ||<1e-4 means convergence and calculation is completed; if it does not converge, update F(i,j) of contact unit and non-contact unit m and P(i,j) n , establish the force balance equation, deformation coordination equation and elastic equation of the contact unit again and solve them to obtain the contact unit force and deformation, repeat ⑤-⑦ until convergence.

[0075] ⑧End the calculation.

[0076] ⑨ Set different pressures to recalculate deformation and seepage. The normal stresses used in this calculation are: 0, 1, 2, 3, 4 and 5 MPa. Figure 5 , Figure 6 and Figure 7 They are the opening, flow rate and pressure distribution under different normal stresses respectively.

[0077] Beneficial effects of this embodiment:

[0078] This embodiment can solve the technical problems of high mesh requirements, long solution time and non-convergence in the prior art. This embodiment performs stress-seepage coupling simulation on a certain crack to obtain the opening distribution diagram, flow velocity diagram and pressure diagram under different normal stresses. This embodiment has fast calculation speed, simple mesh division, fast division speed, and realizes computer automatic processing. This embodiment takes into account the influence of crack contact on seepage, making the calculation more reasonable.

[0079] The above is only a preferred specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A single joint stress-seepage coupling numerical calculation method, characterized in that: The following steps are involved: Obtaining input parameters of single joint stress-seepage coupling, discretizing the input parameters to obtain discretized units, wherein the discretized units include: contact units and non-contact units, and the input parameters include: upper joint coordinates, lower joint coordinates, elevation, initial opening, and contact radius; Based on the contact unit, a contact equation is established, and based on the contact equation, contact stress and deformation are obtained; Based on the deformation, the opening of the non-contact unit is updated to obtain an updated result, and a contact judgment is performed on the updated result. If the contact judgment result is a non-contact unit, the Reynolds equation of the non-contact unit is established, and based on the Reynolds equation, the pressure and flow rate of the non-contact unit are obtained. The convergence is judged based on the pressure. If it does not converge, the contact equation continues to be established until convergence.

2. The single joint stress-seepage coupling numerical calculation method according to claim 1 is characterized in that: The process of discretizing the input parameters includes: The upper joints and lower joints are discretized using rectangular units to obtain discretized units.

3. The single joint stress-seepage coupling numerical calculation method according to claim 1 is characterized in that: The contact equations include: a force balance equation, a deformation coordination equation and an elasticity equation.

4. The single joint stress-seepage coupling numerical calculation method according to claim 1 is characterized in that: Before establishing the contact equation, it also includes: An initialization judgment is performed on the discretization units to determine whether the discretization units are in contact. If so, the discretization units are contact units; if not, the discretization units are non-contact units.

5. The single joint stress-seepage coupling numerical calculation method according to claim 3 is characterized in that: The process of obtaining contact stress and deformation includes: Based on the force balance equation, the Gaussian iteration method is used to solve the force balance equation to obtain the contact stress and deformation of the protrusion.

6. The single joint stress-seepage coupling numerical calculation method according to claim 1, characterized in that: The process of performing contact judgment on the update result includes: If the opening degree of the update result is less than 0, it is a contact unit; if the opening degree of the update result is greater than 0, it is a non-contact unit.

7. The single joint stress-seepage coupling numerical calculation method according to claim 1 is characterized in that: The process of obtaining the pressure and flow rate of the untouched cell includes: Based on the Reynolds equation and the cubic law, the Galerkin method is used to establish an equivalent weak integral form of the Reynolds equation, and a shape function is established using a rectangular unit with bidirectional node interpolation. Based on the shape function and the equivalent weak integral form, the discrete form of each rectangular unit and a set of equations for all unknown node water pressures are obtained. Based on the set of equations, the pressure and flow rate of the uncontacted units are obtained.

8. The single joint stress-seepage coupling numerical calculation method according to claim 1, characterized in that: The process of judging convergence based on the pressure includes: If the opening and pressure of the non-contacting unit satisfy the convergence relationship, it is converged; if the opening and pressure of the non-contacting unit do not satisfy the convergence relationship, it is not converged.

Images