A method for calculating a breaking mechanics model of a hard roof in a deep stope

Abstract

This invention discloses a calculation method for the fracture mechanical model of hard roof in deep mining areas, belonging to the field of mining engineering and strata control technology. Based on the Winkler-Pasternak two-parameter elastic foundation theory, this invention establishes a coal seam-base roof synergistic model, constructs a set of coupled control equations for a multi-layered elastic foundation, and divides the base roof into four mechanical zones: goaf exposure zone, dip support pressure zone, strike support pressure zone, and corner superposition zone. A normal distribution function is used to describe the hump distribution characteristics of the support pressure. Interface continuity and outer boundary convergence conditions are established. The finite difference method is used to discretize the fourth-order partial differential equations into a set of algebraic equations to solve for the deflection field, bending moment field, and shear force field. This invention overcomes the shortcomings of traditional fixed-support boundaries that neglect the deformability of the coal seam, reveals the fracture mechanism of the base roof under elastic foundation boundary conditions, and provides a theoretical basis for the prevention and control of hard roof disasters in deep mining areas.

Description

Technical Field

[0001] This invention belongs to the field of mining engineering and rock strata control technology, specifically relating to a calculation method for a mechanical model of hard roof fracture in deep mining areas. Background Technology

[0002] In deep coal seam mining, the large-scale exposure and sudden fracture of the hard roof are key factors inducing mine pressure disasters. Under the high-stress environment at depth, the hard roof exhibits characteristics of large-span exposure and high energy storage. The elastic energy released by its fracture is transmitted to the working face and surrounding rock of the roadway in the form of stress waves and dynamic loads, which can easily lead to severe mine pressure disasters such as rockbursts and large-area pressure surges. Accurately predicting the location, sequence, and intensity of the hard roof fracture is of great significance for the safe mining of deep resources.

[0003] Existing basic roof mechanical analysis models mainly include rock-beam models and plate models. Among them, the plate model can better reflect the spatial failure characteristics of the roof under bidirectional stress, and is more realistic than the rock-beam model. However, most existing studies simplify the basic roof boundary conditions to fixed or simply supported, assuming that the coal seam has infinite stiffness constraint on the basic roof, and ignoring the deformability and compressibility characteristics of the coal seam and the immediate roof strata. Especially under deep high-stress conditions, the coal body undergoes significant elastoplastic deformation, and its relationship with the basic roof is not an ideal consolidation relationship, but rather there is a certain degree of contact flexibility and deformation coordination effect. In addition, existing elastic ground plate models mostly use single-layer plate models, failing to fully consider the coordinated deformation and interlayer stress transfer between multiple thick and hard rock layers; the bearing pressure distribution mostly adopts uniform or linear assumptions, failing to accurately reflect the hump distribution characteristics of the bearing pressure of the surrounding rock in the goaf; the boundary condition treatment is relatively simplified, lacking fine division of different mechanical regions and strict treatment of interface continuity conditions.

[0004] Therefore, there is an urgent need for a calculation method for the fracture mechanics model of hard roof that can take into account the elastic foundation boundary effect of coal seams, the coordinated deformation of multiple rock layers, the non-uniform support pressure distribution, and refined zoning modeling, in order to solve the technical problem of large deviations between existing theoretical predictions and engineering practice. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method for calculating the mechanical model of fracture of hard roof in deep mining areas.

[0006] To achieve the above objectives, the present invention provides the following technical solution: This application provides a method for calculating the mechanical model of hard roof fracture in deep mining areas, including the following steps: Based on the two-parameter elastic foundation theory, a mechanical model of the synergistic effect between the coal seam and the basic roof is established. A multi-layer rock mass system was constructed using a laminated plate coupling model. Interlayer contact stiffness was introduced to form a set of governing equations describing the bending, shearing and normal coupling of each plate. For the initial fracture condition of the basic roof of the first mining face, the basic roof is divided into multiple mechanical zones, and deflection control differential equations are established for each zone. A continuous function with hump distribution characteristics is used to describe the bearing pressure distribution of the surrounding rock in the goaf, and interface continuity conditions between each zone and outer boundary convergence conditions far from the mining-affected area are established. A finite difference scheme is constructed using the finite difference method to discretize the deflection control differential equations of each partition into a system of algebraic equations. All difference equations are assembled into a large sparse linear algebraic equation system, and the deflection values ​​of each node are obtained by solving the system. Based on the obtained deflection field, the bending moment field is calculated using the differential relationship between bending moment and deflection, and the shear force field is calculated through the bending moment field, thus realizing the numerical solution of the basic top deflection field, bending moment field and shear force field.

[0007] Furthermore, in the aforementioned two-parameter elastic foundation theory, the foundation reaction force p and the plate deflection w and their Laplace operator... 2 w satisfies p=k . wG p 2 w, where k is the Winkler foundation stiffness coefficient, G p This refers to Pasternak shear stiffness.

[0008] Furthermore, in the laminated plate coupling model, the interlayer contact stress p c satisfy: p c =k c (w i- w i+1 )-G c 2 (w i- w i+1 ) In the formula k c G is the interlaminar normal contact stiffness coefficient. c w is the interlaminar shear stiffness coefficient. i w i+1 Let be the deflections of the i-th layer and the (i+1)-th layer, respectively; The governing equations are: A 4 W+B 2 W+C W=Q in 4 For a bitone operator, 2 Let W be the Laplace operator, Q be the deflection vector, A be the bending stiffness matrix, B be the shear stiffness matrix, and C be the normal stiffness matrix.

[0009] Furthermore, the multiple mechanical zones include the goaf exposed zone, the dip support pressure zone, the strike support pressure zone, and the corner superimposed zone; wherein the deflection control differential equation of the goaf exposed zone is a homogeneous equation, and the control equations of the dip support pressure zone, the strike support pressure zone, and the corner superimposed zone all include external load terms determined by the support pressure distribution function.

[0010] Furthermore, the continuous function exhibiting the hump distribution characteristic is a normal distribution function, and its expression is: In the formula, q1 is the load above the goaf, q2 is the far-field load, λ is the stress concentration factor, γ is the unit weight of the overlying strata, h is the coal seam thickness, and a0 is the width of the limit equilibrium zone of the surrounding rock along the strike of the goaf.

[0011] Furthermore, the interface continuity condition includes: at the interface between the exposed area and the support pressure area, deflection continuity and rotation continuity are satisfied; the outer boundary convergence condition is: in the region far from the influence of mining, the deflection is zero and the first partial derivative of the deflection is zero.

[0012] Furthermore, the finite difference method employs a 13-node difference scheme, where the finite difference approximation expression for the biharmonic operator is: In the formula, h is the grid step size, w0 is the deflection of the center node, w1 to w4 are the deflections of adjacent nodes, w5 to w8 are the deflections of diagonal nodes, and w9 to w 12 For the deflection of the remote node.

[0013] Furthermore, the differential relationship between the bending moment field and the deflection field is as follows: In the formula, Mx and My are the bending moments about the x-axis and y-axis, respectively, M xy Let D be the torque, D be the bending stiffness of the thin plate, and μ be Poisson's ratio; the shear field is obtained by calculating the difference between bending moments. , .

[0014] Furthermore, the Winkler foundation stiffness coefficient k ranges from 100 MN / m to 7000 MN / m, and the Pasternak shear stiffness G... p The value ranges from 0 MN / m to 400 MN / m; The elastic modulus E of the basic roof is 30 GPa, the Poisson's ratio μ is 0.25, the width a0 of the limiting equilibrium zone of the surrounding rock along the strike of the goaf and the width b0 of the limiting equilibrium zone along the dip are both 8m, the coal seam thickness h0 is 6m, the thickness of the immediate roof h1 is 4m, and the thickness of the basic roof h2 is 12m.

[0015] Furthermore, by analyzing the relationship between the elastic foundation coefficient k and the Pasternak shear stiffness G... p Impact on the basic top break behavior; Determine the fracture pattern: When k < 1000MN / m, the coal wall area ahead of the main roof fractures before the middle section; when k > 5000MN / m, the middle section fractures first. When G p When the strength is less than 200 MN / m, the fracture sequence is highly uncertain; when G p When the coal face length is greater than 250MN / m, the longer side will inevitably break first, and the distance L ahead of the coal face will be greater. c As k decreases or G p It increases as it grows.

[0016] Compared with the prior art, this application has the following beneficial effects: This invention proposes a calculation method for the mechanical model of hard roof failure in deep mining areas. Based on the two-parameter elastic foundation theory, a collaborative model of the coal seam and the basic roof is established, overcoming the shortcomings of traditional fixed-boundary assumptions that ignore the deformability of the coal seam. By constructing a multi-layer elastic foundation coupling model of laminated plates and introducing interlayer contact stiffness, an accurate description of the collaborative deformation of multiple thick and hard rock strata is achieved. By adopting refined zonal modeling, a normally distributed hump pressure function, and strict interface continuity conditions, the problems of inaccurate description of the surrounding rock pressure distribution in the goaf and simplified boundary treatment in existing models are solved. Combining finite difference numerical solution and bending moment and shear field calculation, the invention quantitatively reveals for the first time the regulation mechanism of elastic foundation stiffness coefficient and shear stiffness on the basic roof failure position offset, failure mode transformation, and advance coal wall distance. This provides a complete theoretical method and calculation basis for the shift from empirical judgment to parameterized accurate prediction of hard roof disaster prevention and control in deep mining areas. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of an elastic foundation mechanical model.

[0018] Figure 2 This is a schematic diagram of a multi-layer elastic foundation coupling model with laminated plates.

[0019] Figure 3 This is a schematic diagram of the initial fracture zone model of the basic top of the elastic foundation.

[0020] Figure 4 The diagram shows the initial fracture mechanical model of the basic roof of the first mining face, where (a) is a plan view, (b) is a cross-sectional view of section I-I, and (c) is a cross-sectional view of section II-II.

[0021] Figure 5 This is a curve of the supporting pressure of the surrounding rock in the goaf.

[0022] Figure 6 The diagram shows the basic top plate partition model under different boundary conditions.

[0023] Figure 7 This is the basic top deflection cloud diagram obtained in Example 1.

[0024] Figure 8 The maximum principal bending moment M1 distribution cloud map obtained in Example 1 is shown.

[0025] Figure 9 The distribution cloud map of the minimum principal bending moment M3 obtained in Example 1.

[0026] Figure 10 The diagram shows the influence of the elastic foundation coefficient k on the principal bending moment. In the diagram, (a) shows the curves of the principal bending moment and the distance to the coal face ahead as a function of k, and (b) shows the curve of the bending moment ratio as a function of k.

[0027] Figure 11 For Pasternak shear stiffness G p The diagram shows the influence of the principal bending moment on the distance from the leading coal wall as a function of G. p The curve (b) shows the change in the moment ratio as a function of G. p Change curve. Detailed Implementation

[0028] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0029] Furthermore, in this invention, an element referred to as fixed to or disposed on another element may be directly disposed on the other element, or there may be an intermediate element. When an element is considered to be connected to another element, it may be directly connected to the other element, or there may be an intermediate element present simultaneously. The terms vertical, horizontal, left, right, and similar expressions used herein are for illustrative purposes only and do not represent the only possible implementation.

[0030] See Figures 1-11 This application provides a calculation method for the fracture mechanical model of hard roof in deep mining areas, including the following steps 1-7: Step 1. Based on the two-parameter elastic foundation theory, establish a mechanical model of the synergistic effect between the coal seam and the basic roof.

[0031] In practical implementation, the thick, hard rock layer overlying the working surface is considered as an elastic thin plate, and its fracture problem is simplified to a small deflection bending problem of an elastic rectangular thin plate; the ratio of the plate thickness h to the short side length l satisfies 1 / 100 ≤ h / l ≤ 1 / 5; let the mid-surface of the thin plate be the xoy plane, the z-axis be perpendicular to the mid-surface and downward, and the deflection of any point on the mid-surface be w(x, y), which displaces only along the z-direction; the displacement components of any point p(x, y, z) within the plate are: The formula for the foundation reaction force in the Winkler-Pasternak model is: Where k is the Winkler foundation stiffness coefficient, with a value ranging from 100 to 7000 MN / m, and G... p The Pasternak shear stiffness has a value ranging from 0 to 400 MN / m. 2 w is the Laplace operator. The governing equation for the bending of a thin plate on a two-parameter elastic foundation is: Where D is the bending stiffness of the thin plate, expressed as D=Eh³ / [12(1-μ 2 E is the basic top elastic modulus, which is 30 GPa, h is the basic top thickness, and μ is Poisson's ratio, which is 0.25.

[0032] Step 2. Construct a multilayer rock mass system with a coupled plate model, introduce interlayer contact stiffness, and form a set of governing equations describing the bending, shearing and normal coupling of each plate. Consider an overburden system consisting of n layers, where each layer is treated as a thin slab on an elastic foundation, such as... Figure 1 As shown. The layers interact with each other through interlayers of weaker lithology or through direct contact; this interaction can be simplified to interlayer contact stress p.c ,like Figure 2 As shown. The interlayer contact stress is modeled using a bilinear spring-shear layer model: Where k c G is the interlaminar normal contact stiffness coefficient. c w is the interlaminar shear stiffness coefficient. i w i +1 represent the deflection of the i-th and (i+1)-th layers, respectively.

[0033] Each slab is placed on a Winkler-Pasternak two-parameter elastic foundation with ground reaction force For the i-th layer slab, its bending control equation must simultaneously consider the foundation reaction force, interlayer contact stress, and external loads: After simplification, the governing equations for the i-th layer are obtained: Define the deflection vector W = [w1, w2, ..., w n ] and load vector Q=[q1,q2,…,q n ] The governing equations of the n-layer plate system are expressed in matrix form as follows: In the formula: A, B, and C are the bending stiffness matrix, shear stiffness matrix, and normal stiffness matrix of the governing equation of the elastic substrate, respectively, and their expressions are: This matrix equation provides a unified description of the bending, shearing, and normal coupling behavior of multilayer overburden systems.

[0034] Step 3. For the initial fracture condition of the basic roof of the first mining face, the basic roof is divided into multiple mechanical zones, and the deflection control differential equations for each zone are established respectively. For the initial failure of the basic roof in the first mining face, the basic roof is divided into four mechanical zones: the goaf exposed zone S1, the dip support pressure zone S2, the strike support pressure zone S3, and the corner superimposed zone S4. Figure 3 As shown. A mechanical model for the initial fracture of the basic roof of the first mining face under elastic boundary conditions is established, as follows: Figure 4 As shown.

[0035] A1A2A3A4 represents the boundary of the goaf before the initial failure of the basic roof. A1A4 has a length of 2a (20m) along the strike of the working face and a length of 2b (60m) along the dip of the working face. D1D2D3D4 represents the boundary of an infinitely elastic plate. h0, h1, and h2 represent the coal seam thickness (6m), the immediate roof thickness (4m), and the basic roof thickness (12m), respectively. a0 and b0 represent the widths of the limiting equilibrium zones of the surrounding rock along the strike and dip of the goaf (both taken as 8m). The deflection control equation for the exposed area of ​​the goaf in region S1 is: Adopting such Figure 5 The normal distribution function shown describes the support pressure curve of the surrounding rock in the goaf: Where q1 is the load above the goaf (0.3MPa), q2 is the far-field load (0.8MPa), λ is the stress concentration factor (value 2.5), γ is the unit weight of the overlying strata, and h is the coal seam thickness.

[0036] The deflection control equation for the inclined support pressure zone in region S2 is: The governing equations for regions S3 and S4 are established using a similar method, differing only in the load surface equations. For the periodic failure condition, the already failed basic roof on the goaf side is given a free boundary condition: For periodic failure conditions, the basic roof that has already failed on the goaf side adopts free boundary conditions, that is, the bending moment and shear force are zero at the cut-off point.

[0037] Step 4. Use a continuous function with hump distribution characteristics to describe the support pressure distribution of the surrounding rock in the goaf, and establish the interface continuity conditions between each zone and the outer boundary convergence conditions far from the mining-affected area. At the boundary between zones, the conditions of geometric compatibility and mechanical continuity must be met, meaning that deflection, rotation, and the resulting bending moment and shear force should remain consistent across the interface. Figure 6 As shown. At the interface between S1 and S2, A1-A4 and A2-A3 maintain continuous deflection and rotation angle: At the interface between S1 and S3, A1-A2 and A3-A4 maintain continuous deflection and rotation angle: At the interface where S2, S3 and S4 intersect, deflection and rotation remain continuous: In areas far from the impact of mining, the deformation converges to zero, i.e., w=0. w / x=0 or w / y=0; Assuming that the basic top deflection and rotation angle are zero when |x|≥x1 or |y|≥y1, the boundary conditions at infinity are transformed into discrete nodal difference equations within a finite range.

[0038] Step 5. Construct a difference scheme using the finite difference method to discretize the deflection control differential equations of each partition into a system of algebraic equations; A square grid with equal step size Δx=Δy=h is used, with 13 nodes including: w0 as the center node, w1~w4 as adjacent nodes, w5~w8 as diagonal nodes, and w9~w... 12 For the distant node. The finite difference approximation of the Laplace operator is: The finite difference approximation of the biharmonic operator is: Example of difference equations for each region: The difference equation for the exposed area S1 in the goaf is: The difference equation for the bearing pressure zone S2 is: Since the impact of mining activities on the rock strata is limited, it is assumed that the basic top deflection and rotation angle are zero when |x|≥x1 or |y|≥y1. The boundary conditions at infinity are transformed into discrete nodal difference equations within a finite range (see Table 1).

[0039] Table 1. Difference Expressions for Boundary Conditions The difference expressions for the boundary conditions are performed according to Table 1: w0=0 and (w1-w3) / (2h)=0 on the boundaries B1-B4 of region S2; w0=0 and (w2-w4) / (2h)=0 on the boundaries C1-C2 of region S3; and similarly on the boundaries of region S4.

[0040] Step 6. Assemble all the difference equations into a large sparse linear algebraic equation system and solve for the deflection values ​​at each node. The difference equations of all internal and boundary nodes are assembled into a large sparse linear algebraic equation system according to the global node numbering: K·w=f, where K is the coefficient matrix (an asymmetric sparse matrix), w is the vector composed of the deflections of all nodes, and f is the load vector. In practice, the matrix assembly is implemented using Matlab software, and the equation system is solved using direct methods (such as the backslash operator "\" or LU decomposition) to obtain the deflection value w(x,y) of each grid node.

[0041] Step 7. Based on the deflection field obtained by the solution, calculate the bending moment field using the differential relationship between bending moment and deflection, and calculate the shear force field through the bending moment field, so as to realize the numerical solution of the basic top deflection field, bending moment field and shear force field.

[0042] After obtaining the deflection at each node, the transformation from deflection to bending moment is established using the differential relationship between the bending moment field and the deflection field: Shear force is obtained from moment difference: , The deflection field, bending moment field, and shear force field distribution of the basic top are obtained.

[0043] This yields the deflection, bending moment, and shear force contour maps for the entire basic roof, allowing us to determine the location of maximum deflection, the extreme values ​​and locations of principal bending moments, the fracture sequence (whether the long-side leading coal face fractures first or the middle section fractures first), and the distance L from the leading coal face. c .

[0044] In this embodiment, the technical solution of this application is adopted, which has the following beneficial effects: (1) Based on the Winkler-Pasternak dual-parameter elastic foundation theory, the shortcomings of the traditional fixed boundary assumption that ignores the deformable characteristics of coal seams are overcome, and an elastic foundation mechanical model considering the synergistic effect of coal seam-basic roof is established, which is more in line with the actual geomechanical characteristics.

[0045] (2) The governing equations of the coupled model of multilayer elastic foundation of laminated plate were derived, and the interlayer normal contact stiffness k was introduced. c and shear stiffness G c The coupled bending problem of an n-layer overburden system is described in a unified matrix form, providing a theoretical framework for the collaborative deformation analysis of multi-layer thick and hard rock strata.

[0046] (3) The basic roof is divided into four mechanical zones: the goaf exposed zone S1, the dip support pressure zone S2, the strike support pressure zone S3, and the corner superposition zone S4. The deflection control differential equations for each zone are established. The hump distribution characteristics of the surrounding rock support pressure in the goaf are described by the normal distribution function. The system of regional interface continuity conditions and outer boundary convergence conditions is established, and the fine division of different mechanical zones is realized.

[0047] (4) A 13-node difference scheme was constructed using the finite difference method, and the fourth-order partial differential equation was discretized into a sparse linear algebraic equation system. The numerical solution of the deflection field, bending moment field and shear force field of the basic top of the elastic foundation was realized. It was revealed that the maximum deflection of the basic top under the boundary conditions of the elastic foundation appears in the center of the goaf and increases significantly with the decrease of the elastic foundation coefficient.

[0048] (5) The elastic foundation stiffness coefficient k and Pasternak shear stiffness G were revealed through parameter sensitivity analysis. p Control mechanism for basic jacking failure behavior: When k increases from 100 to 7000MN / m, the central principal bending moment M z The decrease reached 85%, and the distance from the coal face ahead was L. c The fracture pattern decreased from 20-28m to 1-2m, and the fracture mode changed from edge-first to center-first; G p The impact on the overall fracture mode is limited in the range of 0~200MN / m, but beyond 250~300MN / m, M c A leapfrog growth occurred, L c The depth was increased to 25-30m, providing a theoretical basis for disaster prevention and control of hard roof in deep mining areas.

[0049] In one specific implementation, the two-parameter elastic foundation theory relates the foundation reaction force p to the plate deflection w and its Laplace operator. 2 w satisfies p=k . wG p 2 w, where k is the Winkler foundation stiffness coefficient, G p This refers to Pasternak shear stiffness.

[0050] In practice, after the initial fracture of the basic roof, as the working face continues to advance, the fractured basic roof on one side of the goaf no longer has bearing capacity. At this point, this side boundary is treated as a free boundary. Mathematically, the free boundary condition requires that the bending moment and shear force at the cutoff line be zero: M x =0, M y =0, Q x =0, Q y=0; In the finite difference method, virtual nodes are set outside the boundary, and the difference expression of the free boundary (e.g., M) is used. x =0 gives the relationship of the second derivative) to establish the relationship between the deflection of virtual nodes and the deflection of internal nodes, thereby eliminating virtual nodes and forming the algebraic equation of the boundary nodes; the division and control equation of the remaining zones (goaf exposed zone, front support pressure zone, corner superposition zone) are the same as the initial fracture, but the goaf boundary is changed from fixed support or elastic constraint to free; this implementation method can accurately simulate the roof fracture behavior during the periodic pressure process.

[0051] In one specific embodiment, in the laminated plate coupling model, the interlayer contact stress p c satisfy: p c =k c (w i- w i+1 )-G c 2 (w i- w i+1 ) In the formula k c G is the interlaminar normal contact stiffness coefficient. c w is the interlaminar shear stiffness coefficient. i w i+1 Let be the deflections of the i-th layer and the (i+1)-th layer, respectively; The governing equations are: A 4 W+B 2 W+C W=Q in 4 For a bitone operator, 2 Let W be the Laplace operator, Q be the deflection vector, A be the bending stiffness matrix, B be the shear stiffness matrix, and C be the normal stiffness matrix.

[0052] In practice, the grid step size h is determined based on the computational accuracy and the size of the solution domain; in this example, h = 2m. The positions of the 13 nodes are defined as follows: the center node is the current computation point; w1, w2, w3, and w4 are its four adjacent nodes (distance h) to the left, top, right, and bottom, respectively; w5, w6, w7, and w8 are its four diagonal nodes (distance √2h) to the top left, top right, bottom right, and bottom left, respectively; w9, w 10 w 11 w 12The four distal nodes are designated as left-left, top-up, right-right, and bottom-down (distance 2h). This difference scheme has a truncation error of O(h²) for the fourth-order biharmonic operator, offering higher accuracy compared to the 5-node scheme (using only adjacent nodes), and is particularly suitable for capturing local peak values ​​of bending moment. In the programming implementation, this formula is used to calculate the value for each internal node. 4 w0, and substitute it into the governing equation to obtain the algebraic equation.

[0053] In one specific embodiment, the plurality of mechanical zones include a goaf exposed zone, a dip support pressure zone, a strike support pressure zone, and a corner superimposed zone; wherein the deflection control differential equation of the goaf exposed zone is a homogeneous equation, and the control equations of the dip support pressure zone, the strike support pressure zone, and the corner superimposed zone all include external load terms determined by the support pressure distribution function.

[0054] In practical implementation, for an overburden system with n layers, the elastic modulus E of each layer is first determined. i and thickness h i Calculate the bending stiffness D of each layer of the plate. i =E i h i ³ / [12(1-μ i ²)], forming a diagonal matrix A=diag(D1,D2,…,D n The shear stiffness matrix B is a tridiagonal matrix: the main diagonal elements are the total shear stiffness of the i-th layer (including the foundation shear stiffness G). p,i and interlayer shear stiffness G c,i-1 +G c,i The secondary diagonal element is -G. c,i The normal stiffness matrix C is similar, with the main diagonal elements representing the total normal stiffness of the i-th layer (including the foundation stiffness k). i and the contact stiffness k between upper and lower layers c,i-1 +k c,i The secondary diagonal element is -k. c,i In numerical solutions, this set of matrix equations can be regarded as the governing equations of a generalized elastic substrate, and the same finite difference strategy as for single-layer slabs can be adopted. 4 and ² After discretization, a system of algebraic equations is formed and solved.

[0055] In one specific embodiment, the continuous function exhibiting hump distribution characteristics is a normal distribution function, the expression of which is: In the formula, q1 is the load above the gob area, q2 is the far-field load, λ is the stress concentration coefficient, γ is the unit weight of the overlying strata, h is the coal seam thickness, and a0 is the width of the ultimate equilibrium zone of the surrounding rock along the strike of the gob area.

[0056] In specific implementation, first, determine the load q1 above the gob area, the far-field load q2, the stress concentration coefficient λ, the unit weight γ of the overlying strata, the coal seam thickness h, and the width a0 of the ultimate equilibrium zone according to the mine geological data. Then substitute these parameters into the normal distribution function: q(x)=q2+(λγh - q2)·[(q1 - q2) / (λγh - q2)] (|x|-a 0 )² / a 0 ² .

[0057] When |x| ≤ a0, the exponential part is 0 and the base is 1, and at this time q(x) = q1; when |x| increases, q(x) gradually approaches q2. This function reaches the peak value of λγh at the coal wall front (|x| = a0). It should be noted that it is defined in segments during implementation: for the region where |x| < a0, the abutment pressure has decayed to the gob area load, and in actual calculation, q(x) can be directly set to q1; for the region where |x| ≥ a0, it is calculated according to the above formula. Add this function as the external load term to the right end of the control equations in regions S2, S3, and S4.

[0058] In a specific implementation manner, the interface continuity conditions include: at the interface between the exposed area and the abutment pressure area, the deflection continuity and the rotation angle continuity are satisfied; the outer boundary convergence condition is: in the region far from the mining influence, the deflection is zero and the first-order partial derivative of the deflection is zero.

[0059] During specific implementation: at the interface between S1 and S2 (x = ±b, |y| ≤ a), for each node on the interface, it is required that the deflections are equal: w S1 =w S2 , and at the same time, it is required that the partial derivatives of the deflection along the x direction are equal: w S1 / x= w S2 / x. In the difference format, the rotation angle condition can be expressed as (w S1， right adjacent point - w S1， left adjacent point) / (2h)=(w S2， right adjacent point - w S2， left adjacent point) / (2h); these conditions are used to establish the equations of the nodes on the interface or to eliminate the virtual nodes on both sides of the interface line.

[0060] At the interface between S1 and S3 (y = ±a, |x| ≤ b), similarly, w S1=w S3 and w S1 / y= w S3 / y.

[0061] Convergence condition of the outer boundary: In the region far from the influence of mining (|x|≥x1 or |y|≥y1), it is assumed that the basic top no longer deforms, i.e., w=0, and w / x=0 (for the boundary in the x-direction) or w / y=0 (for the boundary in the y direction). The values ​​of x1 and y1 should be large enough to ensure that the deflection at the boundary has decayed to a negligible level; in this example, 300m is used. In the differential mesh, the nodes on these boundaries are treated as given values ​​(w=0), and the deflection of the virtual nodes outside them is determined by the condition that the first derivative is zero (e.g., for...). w / When x=0, we have w 外侧 =w 内侧 ).

[0062] In one specific implementation, the finite difference method employs a 13-node difference scheme, wherein the finite difference approximation expression of the biharmonic operator is: In the formula, h is the grid step size, w0 is the deflection of the center node, w1 to w4 are the deflections of adjacent nodes, w5 to w8 are the deflections of diagonal nodes, and w9 to w 12 For the deflection of the remote node.

[0063] The bending moment-deflection differential relationship is as follows: In the formula, Mx and My are the bending moments about the x-axis and y-axis, respectively, M xy Where is the torque, D is the bending stiffness of the thin plate, D=Eh2³ / [12(1-μ²)]; μ is Poisson's ratio; the shear field is obtained by calculating the difference of bending moments: , .

[0064] Q y ≈(M y2 -M y4 ) / (2h)+(M xy1 -M xy3 ) / (2h) Here Mx1 M x3 M represents the right and left nodes of the current node, respectively. x Value, M xy2 M xy4 M represents the upper and lower nodes respectively. xy The values ​​are calculated similarly. The bending moment and shear force distributions throughout the entire plate domain can be obtained through the above calculations.

[0065] In one specific embodiment, the Winkler foundation stiffness coefficient k ranges from 100 MN / m to 7000 MN / m, and the Pasternak shear stiffness G... p The value ranges from 0 MN / m to 400 MN / m; The elastic modulus E of the basic roof is 30 GPa, the Poisson's ratio μ is 0.25, the width a0 of the limiting equilibrium zone of the surrounding rock along the strike of the goaf and the width b0 of the limiting equilibrium zone along the dip are both 8m, the coal seam thickness h0 is 6m, the thickness of the immediate roof h1 is 4m, and the thickness of the basic roof h2 is 12m.

[0066] In practice, k is selected between 100 and 7000 MN / m based on the measured coal and rock mechanical parameters in the mine or experience, and G... p Between 0 and 400 MN / m; for typical deep coal seams, k is often taken as 200 to 2000 MN / m, and G... p Take values ​​from 0 to 200 MN / m; in the numerical simulation, select k and G... p Substituting the value into the foundation reaction formula p=kw-G p ²w, which is then reflected in the corresponding terms of the governing equations. By changing k and G p The value of can be used to analyze the influence of different elastic foundation properties on the basic jacking failure behavior (as shown in Examples 2 and 3). This value range covers the vast majority of actual working conditions.

[0067] In practical implementation, the basic elastic modulus E is taken as 30 GPa, and the Poisson's ratio μ is taken as 0.25; the width of the limit equilibrium zone along the strike of the goaf... and the width of the limit equilibrium zone along the dip All are taken as 8m; coal seam thickness Take 6m, directly top thickness Take 4m as the basic top thickness. Take 12m. These parameters are given in Example 1. Users can adjust them according to actual mine geological data, but the above values ​​can be used as typical reference values. Substituting these values ​​into the calculation of bending stiffness D, support pressure function, governing equation coefficients, etc., the specific deflection, bending moment and shear force results can be obtained.

[0068] In one specific implementation, the elastic base coefficient k and Pasternak shear stiffness G are analyzed. p Impact on the basic top break behavior; Determine the fracture pattern: When k < 1000MN / m, the coal wall area ahead of the main roof fractures before the middle section; when k > 5000MN / m, the middle section fractures first. When G p When the strength is less than 200 MN / m, the fracture sequence is highly uncertain; when G p When the coal face length is greater than 250MN / m, the longer side will inevitably break first, and the distance L ahead of the coal face will be greater. c As k decreases or G p It increases as it grows.

[0069] In practical implementation, after completing the modeling and calculation of the fracture mechanics model for hard roofs in deep mining areas, parameter sensitivity analysis is performed. Keeping other parameters constant, the following parameters are varied: k value: (e.g., 100, 500, 1000, 2000, 3000, 5000, 7000MN / m) G p value: (e.g., 0, 25, 50, 100, 150, 200, 250, 300, 350, 400 MN / m); Calculate the principal bending moment distribution for each case. Extract the maximum principal bending moment M in the center of the goaf. z And the minimum principal bending moment M in the area of ​​the coal wall ahead of the main top long side. c (or M) d (Take the larger absolute value) and calculate the ratio M. c / M z Simultaneously, the distance L between the extreme point of the long-side bending moment and the coal face is extracted. c Determine the breaking order based on the ratio: When k < 1000MN / m, M c / M z The value is significantly greater than 1 (approximately 3-5). The bending moment borne by the long side leading coal wall zone is much higher than that in the middle, so the long side breaks first. When k > 5000MN / m, M c / M z When the value is close to or less than 1, the bending moment in the middle becomes dominant, so the middle part fails first. When G p When <200MN / m, M c / M z Fluctuating between 0.7 and 1.1, the bending moments in the middle and the edges are close to competing, and the order of failure is highly uncertain. When G pWhen M >250MN / m, c / M z When the bending moment increases sharply to over 2.0, the bending moment on the longer side becomes significantly dominant, so the longer side will inevitably break first, and at this point, L... c As k decreases or G p Increases as k increases (e.g., L increases when k=100MN / m) c =20~28m, k=7000MN / m when L c =1~2m; G p =0 L c =2~3m, G p =300MN / m when L c =25~30m).

[0070] This judgment rule is based on coal seam conditions (equivalent k and G) in engineering practice. p It provides a quantitative basis for predicting the location and sequence of roof failure, which helps to take targeted pressure relief or support measures.

[0071] Example 1: Mechanical calculation of the initial failure of the basic roof of the first mining face based on the Winkler-Pasternak two-parameter elastic foundation theory.

[0072] Taking a typical deep mine's first mining face as the engineering background, the basic calculation parameters are as follows: basic roof elastic modulus E is 30 GPa, basic roof thickness h2 is 12 m, Poisson's ratio ν is 0.25, face dip half-width b is 60 m, strike advance half-length a is 20 m, load above the goaf q1 is 0.3 MPa, far-field load q2 is 0.8 MPa, coal seam thickness h0 is 6 m, immediate roof thickness h1 is 4 m, elastic foundation stiffness k is 200 MN / m, and foundation shear stiffness G... p The stress concentration factor λ is 2.5, and the widths a0 and b0 of the limit equilibrium zone are both 8m.

[0073] Step 1: Establish the bending control equation for a thin plate on a two-parameter elastic foundation. The bending stiffness of the thin plate is D = Eh. 3 / [12(1-μ 2 )]=30×10 9 ×12 3 / [12(1-0.25 2 )]=4.096×10 11 N·m. Foundation reaction force p=kw(x,y)-G p 2 w=200×10 6 w-50×10 6 2 w; The governing equation is: 4 w / x 4 +2 4 w / ( x 2 y 2 )+ 4 w / y 4 +(50×10 6 / 4.096×10 11 ()( 2 w / x 2 + 2 w / y 2 )=[q(x,y)-200×10 6 w] / 4.096×10 11 ; Step 2: Describe the supporting pressure curve of the surrounding rock in the goaf using a normal distribution function; take γ = 25000 N / m 3 If h = 6m, then: f1(x)=0.8+(2.5×25000×6-0.8×10 6 )×[(0.3-0.8) / (2.5×25000×6-0.8×10 6 )] [(x-8)2 / 64] Step 3: Divide the basic top into four mechanical zones: S1, S2, S3, and S4, and establish the deflection control differential equations for each zone; the control equation for zone S1 (-60≤x≤60m, -20≤y≤20m) is: 4 w1 / x 4 +2 4 w1 / ( x 2 y 2 )+ 4 w1 / y 4 =0.3×10 6 / 4.096×10 11 The governing equations for region S2 (|x|>60m, -20≤y≤20m) are: 4 w / x 4 +2 4 w / ( x 2 y 2 ) + 4 w / y 4 +(50 × 10 6 / 4.096 × 10 11 )( 2 w / x 2 + 2 w / y 2 ) = [f1(|x| - 60) - 200 × 10 6 w] / 4.096 × 10 11 .

[0074] Step 4: Establish the interface continuity conditions; at the interface of x = ±60m (-20 ≤ y ≤ 20m): w S1 / x = w S2 / x, w S1 (±60, y) = w S2 (±60, y) at the interface of y = ±20m (-60 < x < 60m): w S1 / y = w S3 x / y, w S1 (x, ±20) = w S3 (x, ±20).

[0075] Step 5: Set the outer boundary convergence conditions. Assume that when |x| ≥ 300m or |y| ≥ 300m, w = 0 and w / x = 0 or w / y = 0.

[0076] Step 6: Use the finite difference method FDM to construct a 13-node difference format. Select the grid step size h = 2m, generate grid nodes in the solution domain, and discretize the control equations for each region. The difference equation for region S1 is: [w9 + w 10 + w 11 + w 12+2(w5+w6+w7+w8)-8(w1+w2+w3+w4)+20w0](4.096×10 11 / 16)-0.3×10 6 =0; The difference equation for region S2 is: [w9+w 10 +w 11 +w 12 +2(w5+w6+w7+w8)-8(w1+w2+w3+w4)+20w0](4.096×10 11 / 16)+[4w0-8(w1+w2+w3+w4)](50×10 6 / 4)+200×10 6 w0-0.8×10 6 -(2.5×25000×6-0.8×10 6 )×[(0.3-0.8) / (2.5×25000×6-0.8×10 6 )] [(|x|-68)2 / 64] =0.

[0077] Step 7: Assemble all the difference equations into a large sparse linear algebraic equation system K·w=f, and solve it using Matlab software to obtain the deflection values ​​of each node.

[0078] Step 8: Calculate the bending moment field using the bending moment-deflection differential relationship: M x =-4.096×10 11 [2w0-w1-w3+0.25(2w0-w2-w4)] / 4 M y =-4.096×10 11 [2w0-w2-w4+0.25(2w0-w1-w3)] / 4 M xy =-4.096×10 11 ×0.75(-w5+w6+w7-w8) / 16 Step 9: Calculate the shear field: Q x =(M1-M3) / 4, Q y =(M2-M4) / 4 like Figures 7-9 As shown, the contour lines of deflection and bending moment before the initial failure of the basic roof under the elastic foundation boundary condition were obtained. The white dashed rectangle represents the boundary of the goaf.

[0079] The results show that the contour lines of the basic roof deflection are roughly distributed in a concentric ellipse, with the maximum deflection occurring at the center of the goaf, approximately 80 mm. The maximum principal bending moment M1 and the minimum principal bending moment M3 at each node in the middle of the basic roof above the goaf are both positive bending moments. The peak value of the maximum principal bending moment M1 is approximately 280 MN·m, occurring near the center of the goaf. In the area of ​​the coal face ahead of the goaf, both M1 and M3 are negative bending moments. The principal bending moment with the largest absolute value in the long and short sides of the basic roof is the minimum principal bending moment M3, with a peak value of approximately -250 MN·m. Moreover, this extreme value is not located at the coal face boundary, but rather approximately 5 m ahead of the coal face.

[0080] Example 2: Analysis of the Influence of Elastic Foundation Stiffness Coefficient k on the Failure Behavior of the Foundation Top Based on Example 1, keeping other parameters unchanged, the value of the elastic foundation coefficient k was changed, and the main bending moment distribution of the basic roof and the advance coal wall distance were calculated when k=100, 500, 1000, 2000, 3000, 5000, and 7000MN / m respectively. Figure 10 The influence of the Winkler elastic foundation stiffness coefficient k on the main bending moment of the basic roof and the advance coal wall distance is presented.

[0081] The results show that: (1) the principal bending moment in the middle is M z The value decreases monotonically as k increases, from approximately 1050 MN·m when k=100MN / m to approximately 150 MN·m when k=7000MN / m, a decrease of 85%. This is because an increase in k means that the elastic constraint capacity of the coal seam is enhanced, the bearing contribution of the elastic foundation outside the goaf to the basic roof is increased, resulting in a decrease in the effective span length in the middle of the goaf and a reduction in the degree of curvature.

[0082] (2) Principal bending moment M on the long side c and the short side principal bending moment M d The bending moment ratio M increases sharply when k is small, but decreases rapidly as k increases, and tends to level off after k > 2000MN / m. c / M z When k < 1000MN / m, the value reaches as high as 3~5, M d / M z The value is approximately 2, indicating that when the elastic constraint is weak, the bending moment borne by the area of ​​the main roof leading the coal wall is much higher than that in the middle, and it will inevitably break before the middle. When k increases to approximately 2000~3000MN / m, M c / M z When the value drops to approximately 1.0~1.3, approaching the critical transition value, the long side and the middle section may fracture synchronously. When k further increases to above 5000 MN / m, M... c / M z Stable at approximately 1.0~1.1, M d / M zWhen the bending moment drops below 0.6-0.7, the bending moment in the middle becomes dominant, and the basic top tends to break first in the middle.

[0083] (3) Distance L of the coal face ahead c and L d The value of L decreases significantly and monotonically as k increases. When k = 100MN / m, L... c Approximately 20~28m, L d The initial distance is approximately 25-28m; when k=3000MN / m, it decreases to approximately 3-5m; when k=7000MN / m, it approaches 1-2m. The reduction in the distance to the leading coal wall reflects the increased stiffness of the elastic foundation, which makes the bending deformation of the basic roof more concentrated directly above the goaf, and the constraint effect of the coal wall boundary approaches the fixed support boundary condition.

[0084] Example 3: Pasternak shear stiffness G p Analysis of the impact on the basic top fracture behavior Based on Example 1, keeping k=200MN / m and other parameters constant, the Pasternak shear stiffness G is changed. p Calculate G for each of the given values. p Distribution of the main bending moment of the basic roof and the distance ahead of the coal wall when MN / m = 0, 25, 50, 100, 150, 200, 250, 300, 350, 400 MN / m. Figure 11 G is given p The influence of the basic top principal bending moment and failure location.

[0085] The results show that: (1) the principal bending moment in the middle is M z With G p The increase in G shows a slow growth trend, from G p =0, 250 MN·m increases to G p When the density is 300 MN / m, it is approximately 370 MN·m, an increase of about 48%. This indicates that G p The increased horizontal continuity of the foundation allows the elastic foundation near the goaf boundary to more effectively transfer deformation to the surrounding area, thereby increasing the effective span length of the basic roof above the goaf to a certain extent.

[0086] (2) with M z In contrast to the slow growth, M c and M d In G p When G is large, it exhibits a rapid increase. p When M increases from 0 to 200MN / m, c and M d It basically remained within the range of 250~400 MN·m, with little variation; but when G p When it exceeds 250~300MN / m, Mc A sharp increase occurred in G p When the velocity is 300 MN / m, it reaches approximately 1200 MN·m, G p When the value is 400 MN / m, it further increases to approximately 3500 MN·m or more; M d It also exhibits a similar trend of leapfrog growth, but the increase is slightly smaller than M. c This is due to the larger G. p This results in excessive shear continuity in the elastic foundation, causing the foundation top to bear an excessive additional bending moment in the transition zone of the elastic foundation.

[0087] (3) Distance L of the coal face ahead c and L d With G p The increase shows a significant increasing trend, especially in G. p Growth intensifies after reaching >200MN / m; G p =0 L c Approximately 2~3m, G p =300MN / m when L c The shear stiffness increases to approximately 25-30 m; this indicates that the increase in Pasternak shear stiffness expands the range of horizontal deformation transmission of the elastic foundation, and the fracture line is significantly pushed away from the coal face.

[0088] (4) Bending moment ratio M c / M z and M d / M z In G p Within a relatively small range (0~200MN / m), it remains stable between 0.7 and 1.1, close to the critical value of 1.0, indicating that the bending moments at the middle and edges are in competition, and the uncertainty of the failure sequence is relatively large; when G p After exceeding 250MN / m, M c / M z It surged to over 2.0-5.0, M d / M z The temperature also increases simultaneously. At this time, the basic top will inevitably break first at the position of the long side ahead of the coal wall, and the break position is far away from the coal wall.

[0089] (5) Comparison with G p =0 and G p The calculation result of 50MN / m shows that the introduction of the Pasternak term makes M c Slightly increased, L c It increases by approximately 1-2 meters, but its impact on the overall fracture mode is relatively limited. However, when G... p When the value is large, the Pasternak effect becomes significant in influencing the location and sequence of fractures.

[0090] The innovation of this invention lies in the construction of an elastic foundation boundary theory system that considers the deformable characteristics of coal seams. Through the synergistic application of a two-parameter foundation model, multi-layer coupled equations, refined zonal modeling, and a 13-node difference scheme, it reveals for the first time the quantitative control mechanism of elastic constraints on the displacement of the basic roof failure position, the change of failure mode, and the advance coal wall distance. This provides theoretical support for the prevention and control of hard roof disasters in deep mining areas to shift from empirical judgment to parameterized accurate prediction, and has significant academic value and engineering application prospects.

[0091] Advantages compared to existing technologies It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0092] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A calculation method for a mechanical model of hard roof fracture in deep mining areas, characterized in that, Includes the following steps: Based on the two-parameter elastic foundation theory, a mechanical model of the synergistic effect between the coal seam and the basic roof is established. A multi-layer rock mass system was constructed using a laminated plate coupling model. Interlayer contact stiffness was introduced to form a set of governing equations describing the bending, shearing and normal coupling of each plate. For the initial fracture condition of the basic roof of the first mining face, the basic roof is divided into multiple mechanical zones, and deflection control differential equations are established for each zone. A continuous function with hump distribution characteristics is used to describe the bearing pressure distribution of the surrounding rock in the goaf, and interface continuity conditions between each zone and outer boundary convergence conditions far from the mining-affected area are established. A finite difference scheme is constructed using the finite difference method to discretize the deflection control differential equations of each partition into a system of algebraic equations. All difference equations are assembled into a large sparse linear algebraic equation system, and the deflection values ​​of each node are obtained by solving the system. Based on the obtained deflection field, the bending moment field is calculated using the differential relationship between bending moment and deflection, and the shear force field is calculated through the bending moment field, thus realizing the numerical solution of the basic top deflection field, bending moment field and shear force field.

2. The calculation method for the mechanical model of hard roof fracture in deep mining areas according to claim 1, characterized in that, In the aforementioned two-parameter elastic foundation theory, the foundation reaction force p and the plate deflection w and their Laplace operator are related. 2 w satisfies p=k . wG p 2 w, where k is the Winkler foundation stiffness coefficient, G p This refers to Pasternak shear stiffness.

3. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 1, characterized in that, In the laminated plate coupling model, the interlayer contact stress p c satisfy: p c =k c (w i- w i+1 )-G c 2 (w i- w i+1 ) In the formula k c G is the interlaminar normal contact stiffness coefficient. c w is the interlaminar shear stiffness coefficient. i w i+1 Let be the deflections of the i-th layer and the (i+1)-th layer, respectively; The governing equations are: A 4 W+B 2 W+C W=Q in 4 For a bitone operator, 2 Let W be the Laplace operator, Q be the deflection vector, A be the bending stiffness matrix, B be the shear stiffness matrix, and C be the normal stiffness matrix.

4. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 1, characterized in that, The multiple mechanical zones include the goaf exposed zone, the dip support pressure zone, the strike support pressure zone, and the corner superimposed zone; the deflection control differential equation of the goaf exposed zone is a homogeneous equation, and the control equations of the dip support pressure zone, the strike support pressure zone, and the corner superimposed zone all contain external load terms determined by the support pressure distribution function.

5. The calculation method for the mechanical model of hard roof fracture in deep mining areas according to claim 1, characterized in that, The continuous function exhibiting the hump distribution characteristic is a normal distribution function, and its expression is: In the formula, q1 is the load above the goaf, q2 is the far-field load, λ is the stress concentration factor, γ is the unit weight of the overlying strata, h is the coal seam thickness, and a0 is the width of the limit equilibrium zone of the surrounding rock along the strike of the goaf.

6. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 1, characterized in that, The interface continuity conditions include: at the interface between the exposed area and the support pressure area, deflection continuity and rotation continuity are satisfied; the outer boundary convergence condition is: in the area far from the influence of mining, the deflection is zero and the first partial derivative of the deflection is zero.

7. The calculation method for the mechanical model of hard roof fracture in deep mining areas according to claim 1, characterized in that, The finite difference method employs a 13-node difference scheme, where the finite difference approximation expression for the biharmonic operator is: In the formula, h is the grid step size, w0 is the deflection of the center node, w1 to w4 are the deflections of adjacent nodes, w5 to w8 are the deflections of diagonal nodes, and w9 to w 12 For the deflection of the remote node.

8. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 1, characterized in that, The differential relationship between the bending moment field and the deflection field is: In the formula, Mx and My are the bending moments about the x-axis and y-axis, respectively, M xy Let D be the torque, D be the bending stiffness of the thin plate, and μ be Poisson's ratio; the shear field is obtained by calculating the difference between bending moments. , 。 9. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 2, characterized in that, The Winkler foundation stiffness coefficient k ranges from 100 MN / m to 7000 MN / m, and the Pasternak shear stiffness G... p The value ranges from 0 MN / m to 400 MN / m; The elastic modulus E of the basic roof is 30 GPa, the Poisson's ratio μ is 0.25, the width a0 of the limiting equilibrium zone of the surrounding rock along the strike of the goaf and the width b0 of the limiting equilibrium zone along the dip are both 8m, the coal seam thickness h0 is 6m, the thickness of the immediate roof h1 is 4m, and the thickness of the basic roof h2 is 12m.

10. The calculation method for the fracture mechanical model of hard roof in deep mining areas according to claim 1, characterized in that, By analyzing the elastic foundation coefficient k and Pasternak shear stiffness G p Impact on the basic top break behavior; Determine the fracture pattern: When k < 1000MN / m, the coal wall area ahead of the main roof fractures before the middle section; when k > 5000MN / m, the middle section fractures first. When G p When the strength is less than 200 MN / m, the fracture sequence is highly uncertain; when G p When the coal face length is greater than 250MN / m, the longer side will inevitably break first, and the distance L ahead of the coal face will be greater. c As k decreases or G p It increases as it grows.

Images