Method and system for calculating stress intensity factor of elliptic bending crack of quasi-orthogonal anisotropic rock mass

Abstract

The invention discloses a quasi-orthogonal anisotropic rock mass elliptic bending crack stress intensity factor calculation method, which comprises the following steps: acquiring a tunnel surrounding rock crustal stress field, surrounding rock quasi-orthogonal anisotropic parameters and a crack geometrical shape, and constructing a plane unit model containing an elliptic bending crack; constructing an affine transformation and conformal mapping function, and mapping and transforming the bending crack on the z plane to a plane unit circle; deriving an analytical expression of a reset potential on a zeta plane based on an elastic mechanics complex variable function theory; and deriving a calculation formula of the quasi-orthogonal anisotropy rock mass elliptic bending crack stress intensity factor so as to complete calculation of the quasi-orthogonal anisotropy rock mass elliptic bending crack stress intensity factor. The invention also discloses a system for realizing the quasi-orthotropic rock mass elliptic bending crack stress intensity factor calculation method. According to the method, calculation of the quasi-orthotropic rock elliptic bending crack stress intensity factor is achieved, the reliability is higher, the accuracy is better, and the calculation efficiency is higher.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of civil engineering, and particularly relates to a quasi-orthotropic anisotropic rock mass elliptical curved crack stress intensity factor calculation method and system. BACKGROUND

[0002] In the field of geotechnical engineering, such as railway tunnel construction and slope engineering, rock mass with quasi-orthotropic anisotropic characteristics generally has the mechanical property of equal eigenvalues, and the internal structure of the rock mass often develops curved structural joints and fissures. Under the action of the in-situ stress field, the initiation and propagation of such curved cracks may induce major engineering instability accidents and even geological disasters. Therefore, the accurate calculation of the stress field and stress intensity factor of the curved crack in the orthotropic material not only has important engineering practical value for the safety evaluation and stability control of tunnel engineering and slope engineering, but also has significant theoretical significance for in-depth understanding of the mechanical mechanism of the crack propagation path.

[0003] At present, researchers have used various research methods (including energy release rate method, weight function method, perturbation method and singular integral equation method) to systematically study the stress field distribution of the curved crack in the isotropic surrounding rock under the action of the in-situ stress. However, these methods are mainly applicable to two specific cases: one is the micro-radian curved crack with a very small deflection angle, and the other is the short-range curved crack with a very small length of the curved section; therefore, such schemes cannot be directly applied to the calculation of the stress intensity factor of the elliptical curved crack in the quasi-orthotropic anisotropic rock mass. In addition, although some researchers have proposed schemes for theoretically deriving the stress intensity factor of the curved crack with a short kink section in the anisotropic material by using the perturbation method, or have proposed schemes for establishing an analytical model of the short-range curved crack in the anisotropic material based on the energy release rate theory, or have proposed schemes for successfully solving the stress intensity factor of the edge crack and short kink crack in the orthotropic material by introducing a new type of complex function and improving the Stroh form, the calculation process of such schemes is extremely complex, and the reliability and accuracy also have obvious limitations. SUMMARY

[0004] One of the purposes of the present application is to provide a quasi-orthotropic anisotropic rock mass elliptical curved crack stress intensity factor calculation method with high reliability, good accuracy and relative simplicity.

[0005] The second purpose of the present application is to provide a system for implementing the quasi-orthotropic anisotropic rock mass elliptical curved crack stress intensity factor calculation method.

[0006] The quasi-orthotropic anisotropic rock mass elliptical curved crack stress intensity factor calculation method provided by the present application comprises the following steps:

[0007] S1. Obtain the in-situ stress field of the surrounding rock of a tunnel, the quasi-orthotropic anisotropic parameters of the surrounding rock and the geometric shape of the fissure, and construct a planar unit model containing an elliptical curved fissure;

[0008] S2. Construct an affine transformation and a conformal mapping function to map and transform the curved fissure on the z plane to the planar unit circle;

[0009] S3. Based on the complex variable function theory of elastic mechanics, derive the analytic expression of the complex potential on the ζ plane;

[0010] S4. Derive the calculation formula of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved fissure according to the data information obtained in step S3, so as to complete the calculation of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved fissure.

[0011] The obtaining of the in-situ stress field of the surrounding rock of a tunnel, the quasi-orthotropic anisotropic parameters of the surrounding rock and the geometric shape of the fissure and the construction of the planar unit model containing an elliptical curved fissure in step S1 specifically include the following steps:

[0012] Obtain the in-situ stress field of the surrounding rock of a tunnel and wherein is the horizontal far-field stress, is the vertical far-field stress, is the far-field shear stress;

[0013] Construct a quasi-orthotropic anisotropic planar unit model; wherein the orthotropic elastic constants are represented by E x , E y , G xy and v xy , wherein E x is the elastic modulus in the x direction, E y is the elastic modulus in the y direction, G xy is the shear modulus, and v xy is the Poisson's ratio;

[0014] The model inside includes an elliptical curved fissure, and the horizontal distance between the crack tips is 2a, and the angle between the line connecting the elliptical vertex to the crack tip and the horizontal axis is χ.

[0015] The construction of the affine transformation and the conformal mapping function to map and transform the curved fissure on the z plane to the planar unit circle in step S2 specifically includes the following steps:

[0016] Construct the following affine transformation function to transform the curved fissure on the z plane to the curved fissure on the z1 plane:

[0017] z1=x1+iy1=x+iβy

[0018] Where (x1, y1) is the coordinate of the z1 plane; (x, y) is the coordinate of the z plane; i is the imaginary unit; β is a parameter related to the material properties and satisfies

[0019] Construct the following affine transformation function to transform the elliptical bending crack on the z1 plane to the unit circle on the ζ1 plane:

[0020]

[0021] Where ω1(ζ1) is the mapping function of the elliptical bending crack; R1 is the first intermediate variable, and R1 = a1(1-c1 2 ) -1 / 2 / 2, a1 is the length of the arc crack opening corresponding to the minor axis of the elliptical crack, c1 is the second intermediate variable, and c1=sinη1, η1 is the angle between the line from the vertex of the arc crack corresponding to a1 and b1 to the crack tip and the horizontal axis; V1 is the third intermediate variable, and V1=b1(1-c1 2 ) -1 / 2 / 2, b1 is the opening length of the arc corresponding to the major axis of the elliptical crack, and α1 is the fourth intermediate variable and satisfies tanα1=a1tanχ1 / b1, where χ1 is the angle between the line from the vertex to the crack tip in the z1 plane and the horizontal axis.

[0022] The step S3, based on the complex variable function theory of elastic mechanics, deriving the analytical expression of the reset potential on the ζ plane, specifically includes the following steps:

[0023] After affine transformation of the bending crack in the z1 plane, the analytical expression of the reset potential is expressed as

[0024]

[0025] ψ1(z1)=ψ1(ω1(ζ1))≡ψ1(ζ1)

[0026] Plane unit model is only subject to ground stress and The general form of the reset potential is expressed as

[0027]

[0028] ψ1(ζ1)=S2ω1(ζ1)+ψ0(ζ1)

[0029] In the formula is the first complex function to be solved; S1 is the fifth intermediate variable, and ψ0(ζ1) is the second complex function to be solved; S2 is the sixth intermediate variable, and

[0030] and ψ1(ζ1) satisfies the boundary condition, denoted as

[0031]

[0032] where σ is an arbitrary point on the unit circle; is the third complex function; ω1(σ) is the mapping function of the elliptical curved crack; is the conjugate form of ω1'(σ); is the conjugate form of ; f0(σ) is the boundary equation satisfied by the complex function; is the conjugate form of

[0033] The Cauchy integral is performed on both sides of the equation to obtain

[0034]

[0035] where C is the unit circle; is the fourth complex function; ζ1is an arbitrary point on the ζ plane; ω1(σ) is the mapping function; is the conjugate form of is the conjugate form of ψ0(σ);

[0036] Since there exists and where is the complex function to be solved;

[0037] then the solution of is simplified to the solution of the following two integrals:

[0038]

[0039] where I1(ζ1) is the first integral to be solved; I2(ζ1) is the second integral to be solved;

[0040] For I1(ζ1):

[0041] According to the stress conversion formula, the radial stress and shear stress on the surface of the elliptical crack are denoted as

[0042]

[0043] where σ rr is the radial stress on the crack surface; σ rθ is the shear stress on the crack surface; is the horizontal far-field stress; is the vertical far-field stress; ​​is the far field shear stress; γ is the angle between the outer normal vector at an arbitrary point on the elliptical crack and the x1 axis, and χ is the angle measured from the x1 axis to the radial line connecting an arbitrary point on the elliptical crack to the corresponding elliptical geometric center;

[0044] The radial stress and shear stress are recalculated by the following formula:

[0045]

[0046] The angle χ is converted into the following expression

[0047]

[0048] The far field stress is converted into the equivalent surface stress acting on the elliptical curved crack; the expression of the boundary equation is obtained by integrating along the crack surface

[0049]

[0050] Simplifying the above formula, the integral term is written in the following form:

[0051] ∫e 2iγ dz1=it(P1+iP2)-t(P3+iP4)

[0052]

[0053] In the formula, t is the seventh intermediate variable; P1 is the eighth intermediate variable; P2 is the ninth intermediate variable; P3 is the tenth intermediate variable; P4 is the eleventh intermediate variable;

[0054] In which, the calculation formula of P1, P2, P3 and P4 is:

[0055]

[0056]

[0057] Thus, the expression of I1(ζ1) is obtained as

[0058]

[0059] For I2(ζ1):

[0060] is expressed as

[0061]

[0062] In the formula, R1 is the radius of the corresponding circular arc a1; V1 is the radius of the corresponding circular arc b1; T1 is The expression of the first pole in the unit circle; T2 is ​the second pole point in the unit circle; T3 is the first pole point outside the unit circle; T4 is the second pole point outside the unit circle;

[0063] wherein the expression of T1-T4 is

[0064]

[0065] Write the expression of into the Laurent series expansion, expressed as

[0066]

[0067] wherein N0 is the expression of the pole point taking 0 ; N1(T1) is the expression of the pole point taking T1 ; N2(T2) is the expression of the pole point taking T2 ; T1 is the first pole point in the unit circle; T2 is the second pole point in the unit circle; N(σ) is other terms;

[0068] Since and are holomorphic in the unit circle, the specific expression of N0 and N(σ) does not need to be determined, then is simplified as

[0069]

[0070] wherein the expression of N1(T1) and N2(T2) is

[0071]

[0072] Solving the equations, the expression of I2(ζ1) is

[0073]

[0074] Meanwhile, the following conditions are satisfied:

[0075]

[0076] Substitute σ=T1 -1 and σ=T2 -1 into the condition equation, to obtain

[0077]

[0078] Conjugate processing is performed to obtain

[0079]

[0080] Based on the conjugate form, the simplification is obtained

[0081]

[0082] where TN 11 is the twelfth intermediate variable, TN 12 is the thirteenth intermediate variable, BN1is the fourteenth intermediate variable, TN 21 is the fifteenth intermediate variable, TN 22 is the sixteenth intermediate variable, BN2is the seventeenth intermediate variable;

[0083] TN 11 , TN 12 , BN1, TN 21 , TN 22 and BN2expressions are

[0084]

[0085] where N1is the simplified expression of N1(T1); is the conjugate form of N1; N2is the simplified expression of N2(T2); is the conjugate form of N2; T1is the first pole within the unit circle; T2is the second pole within the unit circle; is the conjugate form of T1; is the conjugate form of T2;

[0086] The comprehensive simplification is obtained

[0087]

[0088] where D 11 is the eighteenth intermediate variable; D 12 is the nineteenth intermediate variable; N 11 is the twentieth intermediate variable; N 12 is the twenty-first intermediate variable; D 21 is the twenty-second intermediate variable; D 22 is the twenty-third intermediate variable; N 21 is the twenty-fourth intermediate variable; N 22 is the twenty-fifth intermediate variable; is the conjugate form of N 11 ; is the conjugate form of N 12 ; is the conjugate form of N 21 ; is the conjugate form of N 22 ;

[0089] D 11 、D 12 、N 11 、N 12 、D 21 、D 22 、N 21 and N 22 The expression is

[0090]

[0091] Where i and j are indicator variables, and their values ​​are 1 and 2 respectively;

[0092] Simplified to get and The expression is

[0093]

[0094] In the formula is the conjugated form of BN1; It is the conjugated form of BN2;

[0095] get The expression is

[0096]

[0097] Step S4, based on the data information obtained in step S3, derives a calculation formula for the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass to complete the calculation of the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass, specifically comprising the following steps:

[0098] The stress intensity factor of the elliptical bending crack in the orthotropic element model is obtained by establishing the relationship between the stress intensity factors on the z plane and the z1 plane.

[0099] In the z-plane, the crack stress intensity factor is expressed as

[0100]

[0101] Where K1 is the I-mode stress intensity factor in the z plane; K II is the type II stress intensity factor in the z plane; z is the physical plane; z1 is the analytical plane;

[0102] Rotate the coordinate system so that the tangent line of any crack tip is parallel to the x-axis, and calculate the expression of the stress intensity factor at the crack tip A in the z1 plane as follows:

[0103]

[0104] Where K Ι1K1is the mode I stress intensity factor in the z1plane; ζ ΙΙ1 K2is the mode II stress intensity factor in the z1plane; ζ 1A ζ is the coordinate at the crack tip A; δ1is the angle between the principal axis x and the tangent line at the crack tip; The differential at the point ζ 1A ; ω1”(ζ 1A ) is the second-order differential of ω1(ζ1) at the point ζ 1A ;

[0105] The expression of the stress intensity factor at the crack tip B in the z1plane is

[0106]

[0107] where ζ 1B is the coordinate at the crack tip B; The differential at the point ζ 1B ; ω1”(ζ 1B ) is the second-order differential of ω1(ζ1) at the point ζ 1B ;

[0108] According to the affine transformation, the following stress field mapping relationship is established between the orthotropic z plane and the isotropic z1plane:

[0109]

[0110] The remote stresses in the orthotropic z plane and the isotropic z1plane have the following relationship:

[0111]

[0112]

[0113] In the z1plane, the coordinate system z 01 and the polar coordinate system are constructed, and the coordinate origins of the two coordinate systems are at the crack tip A, and the crack surface coincides with the x 01 axis; in the z 01 coordinate system, the stress components in the range near the crack tip are expressed by the stress intensity factor as

[0114]

[0115] where σ is the stress field component in the z 01 coordinate system; (r1, α1) is the polar coordinate parameter in the z 01 coordinate system;

[0116] In the z1plane, the stress components in the range near the crack tip are expressed as​​

[0117]

[0118] Where δ1 is the distance between x1 axis and x 01 The angle between the axes;

[0119] In the z plane, a new coordinate system z0 and a polar coordinate system (r, α) are constructed; in the z0 plane, the stress component of the set range near the crack tip is expressed as

[0120]

[0121] Where α is the angle between the x-axis and the x0-axis; s1 is the characteristic root of the coordinate system z1, and

[0122] According to the angle between the z plane and the z0 plane, the stress component in the set range near the crack tip in the z plane is expressed as

[0123]

[0124] Where δ is the angle between the x-axis and the x0-axis;

[0125] Combining the relationship between the stress components in the z plane and the z1 plane, the relationship between the stress intensity factors of the cracks in the z plane and the z1 plane is simplified and expressed as

[0126]

[0127] Where K I is the I-mode stress intensity factor in the z plane; K II is the mode II stress intensity factor in the z plane;

[0128] The obtained parameter K I1 , K II1 Substituting into the above expression, we can obtain the calculation result of the crack stress intensity factor in the z plane.

[0129] The application further provides a system for realizing the calculation method of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved crack, comprising a model construction module, a crack mapping module, a complex potential derivation module and a factor calculation module; the model construction module, the crack mapping module, the complex potential derivation module and the factor calculation module are sequentially connected; the model construction module is used for obtaining the ground stress field of the tunnel surrounding rock, the quasi-orthotropic anisotropic parameter of the surrounding rock and the crack geometry, constructing a plane unit model containing an elliptical curved crack, and uploading data to the crack mapping module; the crack mapping module is used for constructing an affine transformation and a conformal mapping function according to the received data information, mapping and transforming the curved crack on the z plane to the plane unit circle, and uploading data to the complex potential derivation module; the complex potential derivation module is used for deriving the analytical expression of the complex potential on the ζ plane based on the complex variable function theory of the elasticity mechanics according to the received data information, and uploading data information to the factor calculation module; the factor calculation module is used for deriving the calculation formula of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved crack according to the received data information, so as to complete the calculation of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved crack.

[0130] The quasi-orthotropic anisotropic rock mass elliptical curved crack stress intensity factor calculation method and system provided by the application realize the calculation of the stress intensity factor of the quasi-orthotropic anisotropic rock mass elliptical curved crack through modeling, derivation, transformation and calculation of the quasi-orthotropic anisotropic rock mass, and have higher reliability, better accuracy and higher calculation efficiency. BRIEF DESCRIPTION OF DRAWINGS

[0131] Figure 1 It is a method flowchart of the method of the application.

[0132] Figure 2 It is a schematic diagram of the orthotropic anisotropic infinite elastic medium model with an arbitrary size elliptical curved crack under the action of remote stress.

[0133] Figure 3 It is a schematic diagram of the conformal transformation of the elliptical curved crack corresponding to the circular crack to the unit circle.

[0134] Figure 4 It is an affine transformation schematic diagram of the method of the application.

[0135] Figure 5 It is a schematic diagram of the stress component near the crack tip A on the z1 plane.

[0136] Figure 6 It is a schematic diagram of the stress component near the crack tip A on the z plane.

[0137] Figure 7The method of the present invention and the ABAQUS finite element method are used to calculate the I-type (K I ) and type II (K II ) Schematic diagram of the change of stress intensity factor with bending angle.

[0138] Figure 8 Schematic diagram of the functional modules of the system of the present invention. DETAILED DESCRIPTION

[0139] like Figure 1 The method flow diagram of the present invention is shown as follows: The method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass disclosed in the present invention comprises the following steps:

[0140] S1. Obtain the in-situ stress field of the tunnel surrounding rock, the quasi-orthogonal anisotropic parameters of the surrounding rock, and the fracture geometry, and construct a plane unit model containing an elliptical curved fracture. This specifically includes the following steps:

[0141] Obtaining the in-situ stress of the tunnel surrounding rock and in is the horizontal far-field stress, is the vertical far-field stress, is the far-field shear stress;

[0142] Build as Figure 2 The quasi-orthogonal anisotropic planar element model shown in Figure 1 is shown in Figure 2; the orthotropic elastic constant is represented by E x 、E y , G xy and v xy Indicates that E x is the elastic modulus in the x direction, E y is the elastic modulus in the y direction, G xy is the shear modulus, v xy is Poisson's ratio;

[0143] The model includes an elliptical curved crack, the horizontal distance between the crack tips is 2a, and the angle between the line from the ellipse vertex to the crack tip and the horizontal axis is χ.

[0144] S2. Construct an affine transformation and conformal mapping function to transform the bending crack on the z plane to the unit circle; specifically, the steps include:

[0145] like Figure 3 As shown, the following affine transformation function is constructed to transform the bending crack on the z plane to the bending crack on the z1 plane:

[0146] z1=x1+iy1=x+iβy

[0147] Where (x1, y1) is the coordinate of the z1 plane; (x, y) is the coordinate of the z plane; i is the imaginary unit; β is a parameter related to the material properties and satisfies

[0148] Construct the following affine transformation function to transform the elliptical bending crack on the z1 plane to the unit circle on the ζ1 plane:

[0149]

[0150] Where ω1(ζ1) is the mapping function of the elliptical bending crack; R1 is the first intermediate variable, and R1 = a1(1-c1 2 ) -1 / 2 / 2, a1 is the arc crack opening length corresponding to the minor axis of the elliptical crack, c1 is the second intermediate variable, and c1 = sinη1, such as Figure 4 As shown, η1 is the angle between the line from the vertex of the arc crack to the crack tip corresponding to a1 and b1 and the horizontal axis; V1 is the third intermediate variable, and V1 = b1 (1-c1 2 ) -1 / 2 / 2, b1 is the opening length of the arc corresponding to the major axis of the elliptical crack, and α1 is the fourth intermediate variable and satisfies tanα1=a1tanχ1 / b1, where χ1 is the angle between the line from the vertex to the crack tip in the z1 plane and the horizontal axis.

[0151] S3. Based on the theory of complex variable functions of elasticity, derive an analytical expression for the reset potential on the ζ plane; specifically, the steps include:

[0152] After affine transformation of the bending crack in the z1 plane, the analytical expression of the reset potential is expressed as

[0153]

[0154] ψ1(z1)=ψ1(ω1(ζ1))≡ψ1(ζ1)

[0155] Plane unit model is only subject to ground stress and The general form of the reset potential is expressed as

[0156]

[0157] ψ1(ζ1)=S2ω1(ζ1)+ψ0(ζ1)

[0158] In the formula is the first complex function to be solved; S1 is the fifth intermediate variable, and ψ0(ζ1) is the second complex function to be solved; S2 is the sixth intermediate variable, and

[0159] and the boundary condition satisfied by ψ1(ζ1) is expressed as

[0160]

[0161] where σ is an arbitrary point on the unit circle; is the third complex variable; ω1(σ) is the mapping function of the elliptical curved crack; is the conjugate form of ω1'(σ); is the conjugate form of is the conjugate form of is the conjugate form of f0(σ); f0(σ) is a complex variable satisfying the boundary equation;

[0162] The Cauchy integral is performed on both sides of the equation to obtain

[0163]

[0164] where C is the unit circle; is the fourth complex variable; ζ1is an arbitrary point on the ζ plane; ω1(σ) is the mapping function; is the conjugate form of is the conjugate form of ψ0(σ);

[0165] Since there exists and where is the complex variable to be solved;

[0166] then the solution of is simplified to the solution of the following two integrals:

[0167]

[0168] where I1(ζ1) is the first integral to be solved; I2(ζ1) is the second integral to be solved;

[0169] For I1(ζ1):

[0170] According to the stress conversion formula, the radial stress and shear stress on the surface of the elliptical crack are expressed as

[0171]

[0172] where σ rr is the radial stress on the crack surface; σ rθ is the shear stress on the crack surface; is the horizontal far-field stress; is the vertical far-field stress; ​​is the far-field shear stress; γ is the angle between the external normal vector at any point of the elliptical crack and the x1 axis, and χ is the angle measured from the x1 axis to the radial line connecting any point on the elliptical crack and the geometric center of the corresponding ellipse;

[0173] The radial stress and shear stress are recalculated using the following formula:

[0174]

[0175] Convert the angle χ into the following expression

[0176]

[0177] The far-field stress is converted into the equivalent surface stress acting on the elliptical bending crack; integrating along the crack surface, the boundary equation is obtained as follows:

[0178]

[0179] Simplify the above formula and write the integral term as follows:

[0180] ∫e 2iγ dz1=it(P1+iP2)-t(P3+iP4)

[0181]

[0182] Wherein t is the seventh intermediate variable; P1 is the eighth intermediate variable; P2 is the ninth intermediate variable; P3 is the tenth intermediate variable; P4 is the eleventh intermediate variable;

[0183] Among them, the calculation formulas for P1, P2, P3 and P4 are:

[0184]

[0185]

[0186] So the expression of I1(ζ1) is obtained as

[0187]

[0188] For I2(ζ1):

[0189] Will Expressed as

[0190]

[0191] Where R1 is the arc radius corresponding to a1; V1 is the arc radius corresponding to b1; T1 is the first pole of the above formula within the unit circle; T2 is the second pole of the above formula within the unit circle; T3 is the first pole of the above formula outside the unit circle; T4 is the second pole of the above formula outside the unit circle;

[0192] Among them, the expressions of T1 to T4 are

[0193]

[0194] Will Written as Laurent series expansion, it is expressed as

[0195]

[0196] Where N0 is the value corresponding to the extreme point 0. Expression; N1(T1) is the corresponding value when the extreme point is T1 Expression; N2(T2) is the corresponding value when the extreme point is T2 Expression; T1 is The first pole in the unit circle; T2 is The second pole inside the unit circle; N(σ) is the other term;

[0197] because and It is holomorphic in the unit circle, so there is no need to determine the specific expressions of N0 and N(σ), then Simplified to

[0198]

[0199] The expressions of N1(T1) and N2(T2) are

[0200] By combining the equations, we can get the expression of I2(ζ1) as

[0201]

[0202] At the same time, the following conditions are met:

[0203]

[0204] Set σ = T1 -1 and σ=T2 -1 Substituting the conditional formula, we get

[0205]

[0206] Perform conjugation treatment to obtain

[0207]

[0208] Based on the conjugate form, simplification gives

[0209]

[0210] where TN 11 is the twelfth intermediate variable, TN 12 is the thirteenth intermediate variable, BN1is the fourteenth intermediate variable, TN 21 is the fifteenth intermediate variable, TN 22 is the sixteenth intermediate variable, BN2is the seventeenth intermediate variable;

[0211] TN 11 , TN 12 , BN1, TN 21 , TN 22 and BN2expressions are

[0212]

[0213] where N1is N1(T1); is the conjugate form of N1; N2is N2(T2); is the conjugate form of N2; T1is the first pole within the unit circle; T2is the second pole within the unit circle; is the conjugate form of T1; is the conjugate form of T2;

[0214] Synthesis simplification gives

[0215]

[0216] where D 11 is the eighteenth intermediate variable; D 12 is the nineteenth intermediate variable; N 11 is the twentieth intermediate variable; N 12 is the twenty-first intermediate variable; D 21 is the twenty-second intermediate variable; D 22 is the twenty-third intermediate variable; N 21 is the twenty-fourth intermediate variable; N 22 is the twenty-fifth intermediate variable; is the conjugate form of N 11 ; is the conjugate form of N 12 ; is the conjugate form of N 21 ; is the conjugate form of N 22 ;

[0217] D11 , D 12 , N 11 , N 12 , D 21 , D 22 , N 21 and N 22 The expression of

[0218]

[0219] Where i and j are index variables, and both take values of 1 and 2;

[0220] Simplifying to obtain and The expression of

[0221]

[0222] Where is the conjugate form of BN1; is the conjugate form of BN2;

[0223] Obtaining The expression of

[0224]

[0225] S4. Deriving the calculation formula of the stress intensity factor of the quasi-orthotropic anisotropic rock body elliptical bending crack to complete the calculation of the stress intensity factor of the quasi-orthotropic anisotropic rock body elliptical bending crack according to the data information obtained in step S3; specifically including the following steps:

[0226] The stress intensity factor of the elliptical bending crack in the orthotropic unit model is obtained by establishing the relationship between the stress intensity factor on the z plane and the z1 plane;

[0227] In the z plane, the crack stress intensity factor is represented as

[0228]

[0229] Where K1 is the I-type stress intensity factor in the z plane; K II is the II-type stress intensity factor in the z plane; z is the physical plane; z1 is the analytical plane;

[0230] The coordinate system is rotated so that the tangent of any tip of the crack is parallel to the x axis, and the expression of the stress intensity factor at the crack tip A in the z1 plane is calculated as follows:

[0231]

[0232] Where K Ι1 is the I-type stress intensity factor in the z1 plane; KΙΙ1 is the type II stress intensity factor in the z1 plane; ζ 1A is the coordinate of the crack tip A; δ1 is the angle between the principal axis x and the tangent line of the crack tip; for At point ζ 1A Differential at ω1" (ζ 1A ) is ω1(ζ1) at point ζ 1A The second-order differential at ;

[0233] The expression of the stress intensity factor at the crack tip B in the z1 plane is:

[0234]

[0235] Where ζ 1B is the coordinate of the crack tip B; for At point ζ 1B Differential at ω1" (ζ 1B ) is ω1(ζ1) at point ζ 1B The second-order differential at ;

[0236] According to the affine transformation, the following stress field mapping relationship is established between the orthotropic z plane and the isotropic z1 plane:

[0237]

[0238] The long-range stresses in the orthotropic z-plane and the isotropic z1-plane have the following relationship:

[0239]

[0240]

[0241] like Figure 5 As shown, the coordinate system z is constructed in the z1 plane 01 and polar coordinate systems, the origin of the coordinates of both is at the crack tip A, and the crack surface is 01 Axes coincide; in z 01 In the coordinate system, the stress component within a set range near the crack tip is expressed using the stress intensity factor as

[0242]

[0243] In the formula For z 01 Stress field components in the coordinate system; (r1, α1) is the z 01 Polar coordinate parameters in the coordinate system;

[0244] In the z1 plane, the stress components of a set range near the crack tip are expressed as

[0245]

[0246] Where δ1 is the distance between x1 axis and x 01 The angle between the axes;

[0247] In the z plane, a new coordinate system z0 and a polar coordinate system (r, α) are constructed; in the z0 plane, the stress component of the set range near the crack tip is expressed as

[0248]

[0249] Where α is the angle between the x-axis and the x0-axis; s1 is the characteristic root of the coordinate system z1, and

[0250] According to the angle between the z plane and the z0 plane, the stress component in the set range near the crack tip in the z plane is expressed as

[0251]

[0252] Where δ is the angle between the x-axis and the x0-axis;

[0253] Combining the relationship between the stress components in the z plane and the z1 plane, the relationship between the stress intensity factors of the cracks in the z plane and the z1 plane is simplified and expressed as

[0254]

[0255] Where K I is the I-mode stress intensity factor in the z plane; K II is the mode II stress intensity factor in the z plane;

[0256] The obtained parameter K I1 , K II1 Substituting into the above expression, we can obtain the calculation result of the crack stress intensity factor in the z plane.

[0257] The method of the present invention is further described below with reference to the following embodiments:

[0258] like Figure 7 The following are the results of quasi-orthotropic materials (E x / E y =0.25, 1, 16) at the tip of the elliptical crack, and the changes of the mode I stress intensity factor and mode II stress intensity factor with the bending angle, and the calculation results of the proposed method are compared with the finite element calculation results of ABAQUS.

[0259] pass Figure 6It can be seen that the magnitude of the stress intensity factor is affected by the bending angle, b / a value, E x / E y The results show that: K I and K II The value gradually increases with the bending angle, and suddenly decreases when α approaches 90°, where K II The increase is less than that of K I , the reduction is greater than K I In addition, with the x / E y Increase, K I The value gradually increases and K II Finally, from the overall point of view, the larger the b / a value, the greater the K I The larger the value, the II The smaller the value.

[0260] like Figure 8 The figure shows a schematic diagram of the functional modules of the system of the present invention: the system disclosed in the present invention for realizing the method for calculating the stress intensity factor of the elliptical bending crack of the quasi-orthogonal anisotropic rock mass comprises a model construction module, a crack mapping module, a reset potential derivation module and a factor calculation module; the model construction module, the crack mapping module, the reset potential derivation module and the factor calculation module are connected in series in sequence; the model construction module is used to obtain the in-situ stress field of the tunnel surrounding rock, the quasi-orthogonal anisotropic parameters of the surrounding rock and the geometric shape of the crack, construct a plane unit model containing the elliptical bending crack, and upload the data to the crack mapping module; the crack mapping module is used to calculate the stress intensity factor of the elliptical bending crack according to the received Data information, construct affine transformation and conformal mapping function, transform the bending crack mapping on the z plane to the plane unit circle, and upload the data to the reset potential derivation module; the reset potential derivation module is used to derive the analytical expression of the reset potential on the ζ plane based on the received data information and the complex variable function theory of elastic mechanics, and upload the data information to the factor calculation module, the factor calculation module is used to derive the calculation formula of the stress intensity factor of the elliptical bending crack of the quasi-orthotropic rock mass according to the received data information, so as to complete the calculation of the stress intensity factor of the elliptical bending crack of the quasi-orthotropic rock mass.

Claims

1. A method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass, comprising the following steps: S1. Obtain the in-situ stress field of the tunnel surrounding rock, the quasi-orthogonal anisotropic parameters of the surrounding rock, and the fracture geometry, and construct a plane element model containing an elliptical curved fracture; S2. Construct an affine transformation and conformal mapping function to transform the bending crack on the z plane to the unit circle; S3. Based on the theory of complex variable functions of elasticity, derive the analytical expression of the reset potential on the ζ plane; S4. Based on the data information obtained in step S3, a calculation formula for the stress intensity factor of the elliptical bending crack in the quasi-orthotropic rock mass is derived to complete the calculation of the stress intensity factor of the elliptical bending crack in the quasi-orthotropic rock mass.

2. The method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass according to claim 1 is characterized in that The step S1 of obtaining the in-situ stress field of the tunnel surrounding rock, the quasi-orthogonal anisotropic parameters of the surrounding rock, and the fracture geometry, and constructing a plane unit model containing an elliptical curved fracture, specifically includes the following steps: Obtaining the in-situ stress of the tunnel surrounding rock and in is the horizontal far-field stress, is the vertical far-field stress, is the far-field shear stress; Construct a quasi-orthogonal anisotropic planar element model; the orthotropic elastic constant is represented by E x 、E y , G xy and v xy Indicates that E x is the elastic modulus in the x direction, E y is the elastic modulus in the y direction, G xy is the shear modulus, v xy is Poisson's ratio; The model includes an elliptical curved crack inside, the horizontal distance between the crack tips is 2a, and the angle between the line from the ellipse vertex to the crack tip and the horizontal axis is χ.

3. The method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass according to claim 2 is characterized in that The construction of the affine transformation and conformal mapping function described in step S2, which transforms the bending crack on the z plane to the planar unit circle, specifically includes the following steps: Construct the following affine transformation function to transform the bending crack on the z plane to the bending crack on the z1 plane: z1=x1+iy1=x+iβy Where (x1, y1) is the coordinate of the z1 plane; (x, y) is the coordinate of the z plane; i is the imaginary unit; β is a parameter related to the material properties and satisfies Construct the following affine transformation function to transform the elliptical bending crack on the z1 plane to the unit circle on the ζ1 plane: Where ω1(ζ1) is the mapping function of the elliptical bending crack; R1 is the first intermediate variable, and R1=a1(1-c1 2 ) -1 / 2 / 2, a1 is the length of the arc crack opening corresponding to the minor axis of the elliptical crack, c1 is the second intermediate variable, and c1=sinη1, η1 is the angle between the line from the vertex of the arc crack corresponding to a1 and b1 to the crack tip and the horizontal axis; V1 is the third intermediate variable, and V1=b1(1-c1 2 ) -12 / 2, b1 is the opening length of the arc corresponding to the major axis of the elliptical crack, and α1 is the fourth intermediate variable and satisfies tanα1=a1tanχ1 / b1, where χ1 is the angle between the line from the vertex to the crack tip in the z1 plane and the horizontal axis.

4. The method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass according to claim 3 is characterized in that The step S3, based on the complex variable function theory of elastic mechanics, deriving the analytical expression of the reset potential on the ζ plane, specifically includes the following steps: After affine transformation of the bending crack in the z1 plane, the analytical expression of the reset potential is expressed as ψ1(z1)=ψ1(ω1(ζ1))≡ψ1(ζ1) Plane unit model is only subject to ground stress and The general form of the reset potential is expressed as In the formula is the first complex function to be solved; S1 is the fifth intermediate variable, and ψ0(ζ1) is the second complex function to be solved; S2 is the sixth intermediate variable, and and ψ1(ζ1) satisfy the boundary conditions, expressed as Where σ is any point on the unit circle; is the third complex variable function; ω1(σ) is the mapping function of the elliptical bending crack; is the conjugated form of ω1′(σ); for The conjugated form of for The conjugate form of ; f0(σ) is the boundary equation satisfied by the complex function; Performing the Cauchy integration on both sides of the equation, we get Where C is the unit circle; is the fourth complex variable function; ζ1 is an arbitrary point on the ζ plane; ω1(σ) is the mapping function; for The conjugated form of is the conjugate form of ψ0(σ); Due to the existence and in is the complex function to be solved; Then the solution will be Simplified to solve the following two integrals: Where I1(ζ1) is the first integral to be calculated; I2(ζ1) is the second integral to be calculated; For I1(ζ1): According to the stress conversion formula, the radial stress and shear stress on the elliptical crack surface are expressed as Where σ rr is the radial stress on the crack surface; σ rθ is the shear stress on the crack surface; is the far-field stress in the horizontal direction; is the far-field stress in the vertical direction; is the far-field shear stress; γ is the angle between the external normal vector at any point of the elliptical crack and the x1 axis, and χ is the angle measured from the x1 axis to the radial line connecting any point on the elliptical crack and the geometric center of the corresponding ellipse; The radial stress and shear stress are recalculated using the following formula: Convert the angle χ into the following expression The far-field stress is converted into the equivalent surface stress acting on the elliptical bending crack; integrating along the crack surface, the boundary equation is obtained as follows: Simplify the above formula and write the integral term as follows: Wherein t is the seventh intermediate variable; P1 is the eighth intermediate variable; P2 is the ninth intermediate variable; P3 is the tenth intermediate variable; P4 is the eleventh intermediate variable; Among them, the calculation formulas for P1, P2, P3 and P4 are: So the expression of I1(ζ1) is obtained as For I2(ζ1): Will Expressed as Where R1 is the arc radius corresponding to a1; V1 is the arc radius corresponding to b1; T1 is The expression of is the first pole in the unit circle; T2 is The expression of is the second pole in the unit circle; T3 is The expression of is the first pole outside the unit circle; T4 is The expression of is at the second pole outside the unit circle; Among them, the expressions of T1 to T4 are Will Written as Laurent series expansion, it is expressed as Where N0 is the value corresponding to the extreme point 0. Expression; N1(T1) is the corresponding value when the extreme point is T1 Expression; N2(T2) is the corresponding value when the extreme point is T2 Expression; T1 is The first pole in the unit circle; T2 is The second pole inside the unit circle; N(σ) is the other term; because and It is holomorphic in the unit circle, so there is no need to determine the specific expressions of N0 and N(σ), then Simplified to The expressions of N1(T1) and N2(T2) are By combining the equations, we can get the expression of I2(ζ1) as At the same time, the following conditions are met: Set σ = T1 -1 and σ=T2 -1 Substituting the conditional formula, we get Perform conjugation treatment to obtain Based on the conjugate form, we can simplify to Where TN 11 is the twelfth intermediate variable, TN 12 is the thirteenth intermediate variable, BN1 is the fourteenth intermediate variable, TN 21 is the fifteenth intermediate variable, TN 22 is the sixteenth intermediate variable, and BN2 is the seventeenth intermediate variable; TN 11 TN 12 、BN1、TN 21 TN 22 and BN2 are expressed as Where N1 is the simplified expression of N1(T1); is the conjugated form of N1; N2 is the simplified expression of N2(T2); is the conjugated form of N2; T1 is The first pole in the unit circle; T2 is The second pole inside the unit circle; It is the conjugated form of T1; It is the conjugated form of T2; Comprehensive simplification gives Where D 11 is the eighteenth intermediate variable; D 12 is the nineteenth intermediate variable; N 11 is the twentieth intermediate variable; N 12 is the twenty-first intermediate variable; D 21 is the twenty-second intermediate variable; D 22 is the twenty-third intermediate variable; N 21 is the twenty-fourth intermediate variable; N 22 is the twenty-fifth intermediate variable; N 11 The conjugated form of N 12 The conjugated form of N 21 The conjugated form of N 22 The conjugated form of D 11 、D 12 、N 11 、N 12 、D 21 、D 22 、N 21 and N 22 The expression is Where i and j are indicator variables, and their values ​​are 1 and 2 respectively; Simplified to get and The expression is In the formula is the conjugated form of BN1; It is the conjugated form of BN2; get The expression is 5. The method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass according to claim 4 is characterized in that Step S4, based on the data information obtained in step S3, derives a calculation formula for the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass to complete the calculation of the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass, specifically comprising the following steps: The stress intensity factor of the elliptical bending crack in the orthotropic element model is obtained by establishing the relationship between the stress intensity factors on the z plane and the z1 plane. In the z-plane, the crack stress intensity factor is expressed as Where K1 is the I-mode stress intensity factor in the z plane; K II is the type II stress intensity factor in the z plane; z is the physical plane; z1 is the analytical plane; Rotate the coordinate system so that the tangent line of any crack tip is parallel to the x-axis, and calculate the expression of the stress intensity factor at the crack tip A in the z1 plane as follows: Where K Ι1 is the I-mode stress intensity factor in the z1 plane; K ΙΙ1 is the type II stress intensity factor in the z1 plane; ζ 1A is the coordinate of the crack tip A; δ1 is the angle between the principal axis x and the tangent line of the crack tip; for At point ζ 1A Differential at ω1" (ζ 1A ) is ω1(ζ1) at point ζ 1A The second-order differential at ; The expression of the stress intensity factor at the crack tip B in the z1 plane is: Where ζ 1B is the coordinate of the crack tip B; for At point ζ 1B Differential at ω1" (ζ 1B ) is ω1(ζ1) at point ζ 1B The second-order differential at ; According to the affine transformation, the following stress field mapping relationship is established between the orthotropic z plane and the isotropic z1 plane: The remote stresses in the orthotropic z-plane and the isotropic z1-plane have the following relationship: Construct coordinate system z in z1 plane 01 and polar coordinate systems, the origin of the coordinates of both is at the crack tip A, and the crack surface is 01 Axes coincide; in z 01 In the coordinate system, the stress component within a set range near the crack tip is expressed using the stress intensity factor as In the formula For z 01 Stress field components in the coordinate system; (r1, α1) is the z 01 Polar coordinate parameters in the coordinate system; In the z1 plane, the stress components of a set range near the crack tip are expressed as Where δ1 is the distance between x1 axis and x 01 The angle between the axes; In the z plane, a new coordinate system z0 and a polar coordinate system (r, α) are constructed; in the z0 plane, the stress component of the set range near the crack tip is expressed as Where α is the angle between the x-axis and the x0-axis; s1 is the characteristic root of the coordinate system z1, and According to the angle between the z plane and the z0 plane, the stress component in the set range near the crack tip in the z plane is expressed as Where δ is the angle between the x-axis and the x0-axis; Combining the relationship between the stress components in the z plane and the z1 plane, the relationship between the stress intensity factors of the cracks in the z plane and the z1 plane is simplified and expressed as Where K I is the I-mode stress intensity factor in the z plane; K II is the mode II stress intensity factor in the z plane; The obtained parameter K I1 , K II1 Substituting into the above expression, we can obtain the calculation result of the crack stress intensity factor in the z plane.

6. A system for implementing the method for calculating the stress intensity factor of an elliptical bending crack in a quasi-orthotropic rock mass according to any one of claims 1 to 5, characterized in that It includes a model construction module, a crack mapping module, a reset potential derivation module, and a factor calculation module; the model construction module, the crack mapping module, the reset potential derivation module, and the factor calculation module are connected in series in sequence; the model construction module is used to obtain the ground stress field of the tunnel surrounding rock, the quasi-orthogonal anisotropic parameters of the surrounding rock, and the crack geometry, construct a plane unit model containing elliptical curved cracks, and upload the data to the crack mapping module; The crack mapping module is used to construct an affine transformation and conformal mapping function based on the received data information, transform the bending crack mapping on the z plane to the plane unit circle, and upload the data to the reset potential derivation module; The reset potential derivation module is used to derive the analytical expression of the reset potential on the ζ plane based on the received data information and the theory of complex variable functions of elastic mechanics, and upload the data information to the factor calculation module. The factor calculation module is used to derive the calculation formula of the stress intensity factor of the elliptical bending crack of the quasi-orthotropic rock mass based on the received data information, so as to complete the calculation of the stress intensity factor of the elliptical bending crack of the quasi-orthotropic rock mass.

Images