CN115906704A - Single-joint stress-seepage coupling numerical calculation method - Google Patents

Abstract

The invention discloses a method for calculating a single-joint stress-seepage coupling numerical value, which belongs to the field of fluid-solid coupling numerical value calculation and comprises the following steps: acquiring input parameters of single-joint stress-seepage coupling, and performing discretization processing on the input parameters to obtain a discretization unit, wherein the discretization unit comprises: a contact unit and a non-contact unit; establishing a contact equation based on the contact unit, and obtaining contact stress and deformation based on the contact equation; updating the opening degree of the untouched unit based on the deformation to obtain an updated result, performing contact judgment on the updated result, establishing a Reynolds equation of the untouched unit based on the judged result, obtaining the pressure and the flow velocity of the untouched unit based on the Reynolds equation, judging the convergence based on the pressure, and if the untouched unit does not converge, continuously establishing a contact equation until the contact equation converges. The invention can solve the technical problems of high requirement on grids and long solving time, has high calculation speed and simple grid division, and realizes the automatic processing of a computer.

Description

A numerical method for stress-seepage coupling in single joints

Technical Field

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

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.

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, a long solution time, and often encounters non-convergence problems. Therefore, CFD simulation is more troublesome to apply.

Summary of the invention

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.

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

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;

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. 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.

Preferably, the process of discretizing the input parameters includes:

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.

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

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

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.

Preferably, the process of obtaining the 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.

Preferably, 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.

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

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.

Preferably, 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.

The technical effects of the present invention are:

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 the pressure and flow rate of the non-contact unit are obtained based on the Reynolds equation, and convergence is judged based on the pressure. If it does not converge, the contact equation continues to be established until convergence.

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

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:

FIG1 is a flow chart of a calculation method in an embodiment of the present invention;

FIG2 is a schematic diagram of a contact area and a non-contact area of a discrete crack unit in an embodiment of the present invention;

FIG3 is a schematic diagram of a crack obtained by splitting a rock block in an embodiment of the present invention;

FIG4 is a schematic diagram of the measurement of crack morphology, elevation and opening in an embodiment of the present invention;

FIG5 is a schematic diagram of the opening distribution under no stress in an embodiment of the present invention;

FIG6 is a flow velocity diagram under different stresses in an embodiment of the present invention;

FIG. 7 is a diagram showing pressure distribution under different stresses in an embodiment of the present invention.

DETAILED DESCRIPTION

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.

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.

Embodiment 1

As shown in FIG1 , this embodiment provides a single joint stress-seepage coupling numerical calculation method, including:

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

S2. Discretization: The upper and lower joints are discretized using rectangular units, and each unit is denoted as (i, j), as shown in Figure 2.

S3. Initialization judgment: judge whether each unit is in contact, the contact unit is recorded as (i,j) _m , and the non-contact unit is recorded as (i,j) _n

S4. Establish the force balance equation, deformation coordination equation and elastic equation of the contact area (i, j) _m and solve them to obtain the magnitude of the contact unit force and displacement.

The force balance equation is as follows:

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.

The Gaussian iteration method is used to solve the force balance equations to obtain the stress distribution F(i,j) _m and deformation Δu of the protrusion.

S5. Based on the deformation value Δu obtained in the previous step, the opening b(i, j) _n of the non-contact area is updated using the opening calculation formula.

The calculation formula of opening is as follows:

b(i,j) _n =b(i,j) _n -Δu _k (2)

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 , and the opening is 0; if b(i,j)n> 0, the unit is not in contact, and the non-contact unit is recorded as (i,j) _n

S7. Establish the overall conduction equation of the discretization of the Reynolds equation of the non-contact unit (i, j) _n and solve it to obtain the pressure and flow velocity in the non-contact area.

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

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

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

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

^{KePe ＝ Fe} (7)

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:

KP＝F (11)

常数项

采用预处理共轭梯度法对式(11)进行求解，得到未接触区域单元上的水压力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:

S8. Solve for the pressure value of the uncontacted unit P(i,j) _m

S9. Determine whether it converges: if ||b _k -b _k+1 ||<1e-4, ||P _k -P _k+1 ||<1e-4, then it converges and the calculation is completed, where b _k and b _k+1 represent the opening sizes of the k and k+1 iteration steps of the non-contact unit, respectively, and P _k and P _k+1 represent the pressure sizes of the k and k+1 iteration steps of the non-contact unit, respectively;

If it does not converge, update F(i,j) _m and P(i,j) _n of the contact unit and the non-contact unit, 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.

S10. End calculation.

The specific implementation cases are as follows:

① First, the 150mm×150mm×150mm granite specimen was Brazilian split to obtain fresh and complete cracks, as shown in Figure 3;

②Measurement of crack morphology data, initial opening and contact radius. First, stick marking points on the four sides of the upper and lower crack specimens (as shown in Figure 4), and use a three-dimensional optical scanning system to scan the upper and lower crack rough surfaces and the contours under the "closure" to obtain point cloud data; then import it into the Geomaic studio software, use the "stitching" function of the software to match the upper and lower crack marking points, and "simulate" the crack closure state (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 from each point to obtain the contact radius value.

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

④ Initialization judgment: judge whether each unit is in contact, the contact unit is recorded as (i,j) _m , and the non-contact unit is recorded as (i,j) _n

⑤ Establish and solve the force balance equation, deformation coordination equation and elastic equation of the contact area (i, j) _m to obtain the magnitude of the contact unit force and displacement. At the same time, update the new opening value and determine whether there is contact.

⑥ Establish the overall conduction equation of the discretization of the Reynolds equation of the non-contact unit (i, j) _n and solve it to obtain the pressure and flow velocity of the non-contact area. And solve the pressure value of the non-contact unit.

⑦ Determine whether it converges: If ||b _k -b _k+1 ||<1e-4, ||P _k -P _k+1 ||<1e-4, it converges and the calculation is completed; if it does not converge, update F(i,j) _m and P(i,j) _n of the contact unit and the non-contact unit, 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.

⑧End the calculation.

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

Beneficial effects of this embodiment:

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.

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 (8)

1. A single joint stress-seepage coupling numerical calculation method, characterized in that it comprises the following steps: 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; 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. 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. 2. The single-joint stress-seepage coupling numerical calculation method according to claim 1, wherein the process of discretizing the input parameters comprises: 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. 3. The single-joint stress-seepage coupling numerical calculation method according to claim 1 is characterized in that the contact equation includes: a force balance equation, a deformation coordination equation and an elastic equation. 4. The single joint stress-seepage coupling numerical calculation method according to claim 1, 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 the contact stress and deformation comprises: 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 updated result comprises: 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, characterized in that the process of obtaining the pressure and flow velocity of the uncontact unit comprises: 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 comprises: 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