CN108763725A - 基于采空区压实效应的工作面开采全过程数值模拟方法 - Google Patents
基于采空区压实效应的工作面开采全过程数值模拟方法 Download PDFInfo
- Publication number
- CN108763725A CN108763725A CN201810505702.3A CN201810505702A CN108763725A CN 108763725 A CN108763725 A CN 108763725A CN 201810505702 A CN201810505702 A CN 201810505702A CN 108763725 A CN108763725 A CN 108763725A
- Authority
- CN
- China
- Prior art keywords
- working face
- model
- caving
- goaf
- zone
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 49
- 238000005065 mining Methods 0.000 title claims abstract description 27
- 239000011435 rock Substances 0.000 claims abstract description 56
- 238000004088 simulation Methods 0.000 claims abstract description 20
- 238000004364 calculation method Methods 0.000 claims abstract description 18
- 239000000463 material Substances 0.000 claims abstract description 12
- 230000000737 periodic effect Effects 0.000 claims abstract description 12
- 230000006378 damage Effects 0.000 claims abstract description 9
- 239000003245 coal Substances 0.000 claims description 38
- 238000009412 basement excavation Methods 0.000 claims description 25
- 238000006073 displacement reaction Methods 0.000 claims description 16
- 238000010008 shearing Methods 0.000 claims description 9
- 238000010835 comparative analysis Methods 0.000 claims description 4
- 238000012669 compression test Methods 0.000 claims description 4
- 230000006835 compression Effects 0.000 claims description 3
- 238000007906 compression Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 238000000205 computational method Methods 0.000 claims description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 2
- 238000010276 construction Methods 0.000 abstract description 2
- 238000012821 model calculation Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 2
- 238000005452 bending Methods 0.000 description 2
- 238000003825 pressing Methods 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 238000005303 weighing Methods 0.000 description 2
- 206010044565 Tremor Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- FFBHFFJDDLITSX-UHFFFAOYSA-N benzyl N-[2-hydroxy-4-(3-oxomorpholin-4-yl)phenyl]carbamate Chemical compound OC1=C(NC(=O)OCC2=CC=CC=C2)C=CC(=C1)N1CCOCC1=O FFBHFFJDDLITSX-UHFFFAOYSA-N 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000977 initiatory effect Effects 0.000 description 1
- 229910052742 iron Inorganic materials 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000010206 sensitivity analysis Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/06—Power analysis or power optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Business, Economics & Management (AREA)
- General Physics & Mathematics (AREA)
- Tourism & Hospitality (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Strategic Management (AREA)
- Economics (AREA)
- General Health & Medical Sciences (AREA)
- Human Resources & Organizations (AREA)
- Marketing (AREA)
- Primary Health Care (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Business, Economics & Management (AREA)
- Animal Husbandry (AREA)
- Mining & Mineral Resources (AREA)
- Agronomy & Crop Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Marine Sciences & Fisheries (AREA)
Abstract
本发明公开一种基于采空区压实效应的工作面开采全过程数值模拟方法，包括如下步骤：(1)收集待评估的工作面的相关地质信息；(2)建立数值计算模型；(3)确定工作面初次来压(lp)和周期来压步距(ls)；(4)确定工作面上方垮落带和裂隙带的范围，并将垮落带单元体定位为空模型；(5)判断采空区是否已被完全充填；(6)将垮落带的单元模型由空模型改为双屈服模型，同时弱化裂隙带范围内岩体的材料属性；(7)重复执行上述过程，直至工作面开采结束。本发明数值模拟方法能很好的模拟长壁开采过程，揭示矿压显现规律，为灾害防控提供理论依据和施工指导。
Description
技术领域
本发明涉及一种工作面开采全过程数值模拟方法，具体涉及一种适用于煤矿井下长壁工作面开采过程中影响矿压显现的基于采空区压实效应的数值模拟方法。
背景技术
随着我国煤炭开采深度与难度的不断增大，开采引发的冲击矿压、煤与瓦斯突出和强矿压显现等工程灾害日渐严峻，严重制约了矿井的安全高效生产。因此，研究采场矿压显现规律，为灾害防控提供理论依据和现场指导具有重要的理论和现实意义。
研究表明，冲击矿压的发生主要与工作面采动过程中煤岩体所承载的静载荷和动载荷有关。采掘空间周围煤岩体中的静载荷主要由围岩自重应力、构造应力和采动后覆岩结构演化造成的采动应力等组成，动载扰动主要包括采掘空间顶板的来压、上覆岩层运动、工作面煤体破坏以及断层滑移等导致的采掘煤岩矿震活动。
综上所述，研究采动后覆岩结构的演化特征及其影响范围对于工作面动静载分布有重要意义。而针对采矿等大尺寸、复杂工程问题的研究，数值模拟以其直观、高效、便于反演现场复杂地质条件和进行大量参数敏感性分析的特点与优势，在矿压研究中得到广泛应用。然而，目前的数值模型计算所得结果与现实存有较大差异。开采过程中顶板破裂、垮落和采空区充填、承压等反复卸压、支撑作用对应力分布影响显著。故进行数值软件的二次开发对采场矿压显现规律的研究是非常必要的。
发明内容
本发明的目的在于克服现有技术中的不足，提供一种基于采空区压实效应的工作面开采全过程数值模拟方法，能够给出的指标物理意义明确，可操作性强，且能较好的再现长壁开采过程，有助于揭示矿压显现规律。
为解决现有技术问题，本发明公开了一种基于采空区压实效应的工作面开采全过程数值模拟方法，包括如下步骤：
(1)收集待评估的工作面的相关地质信息；
(2)根据所述地质信息建立数值计算模型；
(3)根据开挖步距和顶板下沉量确定工作面初次来压步距和周期来压步距；
(4)根据计算过程中单元体塑性应变的变化确定工作面上方垮落带和裂隙带的范围，并将垮落带单元体定位为空模型；
(5)根据垮落带单元体的碎胀系数判断采空区是否已被完全充填；
(6)将垮落带的单元模型由空模型改为双屈服模型，同时弱化裂隙带范围内岩体的材料属性；
(7)重复执行上述过程，直至工作面开采结束。
作为优选方案，
步骤(1)中，相关地质信息包括：岩层的岩性、煤岩体的物理力学属性、工作面埋深、工作面长度以及原岩应力场；
煤岩体的物理力学属性包括煤岩体厚度、弹性模量、泊松比、密度、内摩擦角、内聚力和单轴抗拉强度。
作为优选方案，
步骤(2)中，建立数值计算模型的方法包括如下步骤：
建立模型网格，模块参数设置，边界条件与初始条件设置，模型初始平衡。
作为优选方案，
步骤(3)中，确定工作面初次来压和周期来压步距的方法包括如下步骤：
(31)根据数值模型网格的大小确定开挖步距，设定开挖步距为一个网格；
(32)每运行100步，进行最大不平衡力收敛确定，若最大平衡比率小于1e-5，说明模型受力平衡，可进行进一步的开挖计算；若最大平衡比率大于1e-5，说明模型受力不平衡，此时进行下一步(33)的顶板下沉量判断；
(33)每进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移大于等于开挖煤层的厚度，则说明垮落充分，此时选取该开挖步距内任一点的开挖值为初次来压步距，若最大垂直位移小于开挖煤层的厚度，则继续开挖；
(34)以垮落充分的开挖步距为起点，每继续进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移大于等于开挖煤层的厚度，则选取当前开挖步距内任一点的开挖值为周期来压步距，若最大垂直位移小于开挖煤层的厚度，则继续开挖。
作为优选方案，
步骤(4)中，确定工作面上方垮落带和裂隙带的范围的方法包括如下方法：
基于塑性理论，最大主应变ε1 t表达式如下：
ε1 t＝ε1 e+ε1 p；
其中，ε1 e为最大主弹性应变，ε1 p为最大主塑性应变，εc为顶板岩层单元体垮落前的最大拉应变；若单元体的塑性应变值ε1 t＝0时，该单元体属于弯曲下沉带，若0<ε1 t<εc时，该单元体属于裂隙带，若ε1 t≥εc时，该单元体属于垮落带；其中，εc根据实测的地表下沉量与数值模拟结果对比分析确定。
作为优选方案，
步骤(5)中，垮落带单元体的碎胀系数的表达式为：
其中，hc为煤层厚度，单位为m；hcr为垮落带岩层的厚度，单位为m；
采用碎胀系数来确定采空区是否被垮落带岩体充填满，当BF≥1.43时，说明采空区未被充填满；当BF<1.43时，说明采空区已经被充填满。
作为优选方案，
步骤(6)中，材料属性的计算方法包括如下内容：
数值模拟计算中垮落带充填体的相关力学参数的计算主要采用下式：
其中，E0为初始切线模量，MPa；σc为垮落带岩体的单轴抗压强度，MPa；b为垮落带岩体的碎胀系数；
建立一个尺寸为1m×1m×1m的单元体进行单轴压缩试验，在单元体上表面施加一固定速度并固定边界；通过不断调整密度、体积模量、剪切模量等物理参量使其与Salamon理论计算曲线吻合，最终得到符合数值模拟要求垮落带岩体的的力学参数；
Salamon表达式为：
其中，σ为当受到侧向约束时加载在松散岩体上的压应力；E0为松散岩体的初始切向应力；εv为体积应变；为最大体积应变；其中，表达式为：
其中，BF为垮落带单元体的碎胀系数；
通过FLAC3D软件将空模型行改为双屈服模型，然后利用上述步骤计算所得的相应力学参数对双屈服模型进行重新赋值；
将裂隙带岩体可等效为弹性各向同性材料，然后降低其弹性模量和泊松比。
作为优选方案，
步骤(6)中，一般裂隙带岩体的弹性模量为500MPa，泊松比为0.3，相应的单元体体积模量和剪切模量分别为133.36MPa和61.53MPa。
本发明具有的有益效果：通过对开采全过程中煤岩体损伤软化、顶板垮落、采空区充填及承压等物理过程的分析，采用损伤软化模型描述煤岩体损伤软化性质，通过岩体的拉应变判断顶板垮落，利用双屈服模型模拟采空区承压特性并相应的弱化裂隙带的相关力学属性，以此实现对开采全过程中影响矿压显现规律的主要过程的模拟。该发明数值模拟方法能很好的模拟长壁开采过程，揭示矿压显现规律，为灾害防控提供理论依据和施工指导。
附图说明
图1为本发明的流程图；
图2为本发明中基于3108工作面的数值计算模型图；
图3为本发明中的Salamon理论模型和数值模拟模型的比较图；
图4、图5为本发明中的工作面1开挖过程中覆岩垮落带和裂隙带演化特征图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
如图1所示，一种基于采空区压实效应的工作面开采全过程数值模拟方法，包括以下步骤：
(1)收集待评估的工作面的相关地质信息。
所述的相关地质信息包括：岩层的岩性、煤岩体的物理力学属性、工作面埋深、工作面长度以及原岩应力场(水平应力与垂直应力的比值)；所述煤岩体的物理力学属性包括煤岩体厚度、弹性模量、泊松比、密度、内摩擦角、内聚力和单轴抗拉强度。
(2)根据步骤一中收集到的资料建立数值计算模型。
该步骤包括：在FLAC3D软件中建立模型网格，模块参数设置，边界条件与初始条件设置，模型初始平衡。
(3)根据开挖步距和顶板下沉量(dz max)确定工作面初次来压(lp)和周期来压步距(ls)。
其具体过程包括以下步骤：
(31)根据数值模型网格的大小确定开挖步距，设定开挖步距为一个网格；
(32)每运行100步，进行最大不平衡力收敛确定，若最大平衡比率小于1e-5，说明模型受力平衡，可进行进一步的开挖计算；若最大平衡比率大于1e-5，说明模型受力不平衡，此时进行下一步(33)的顶板下沉量判断；
(33)每进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移(即下沉量)大于等于开挖煤层的厚度，则说明垮落充分，此时选取该开挖步距内任一点的开挖值为初次来压步距(lp)，若最大垂直位移小于开挖煤层的厚度，则继续开挖；
(34)以垮落充分的开挖步距为起点，每继续进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移大于等于开挖煤层的厚度，则选取当前开挖步距内任一点的开挖值为周期来压步距(ls)，若最大垂直位移小于开挖煤层的厚度，则继续开挖。
(4)根据计算过程中单元体塑性应变(ε1 t)的变化确定工作面上方垮落带和裂隙带的范围，并将垮落带单元体定位为空模型。
其具体过程包括以下内容：
基于塑性理论，最大主应变ε1 t表达式如下：
ε1 t＝ε1 e+ε1 p；
其中，ε1 e为最大主弹性应变，ε1 p为最大主塑性应变。
由于垮落带岩体通常都是拉伸破坏，因此垮落带内塑性应变值比弹性应变值大的多，即ε1 p>>ε1 e。因此，可认为最大主应变约等于最大塑性应变，即ε1 t≈ε1 p。
数值模型中将顶板岩层单元体垮落前的最大拉应变定义为临界垮落带应变值(εc)，用此临界值确定垮落带的范围。若单元体的塑性应变值ε1 t＝0时，该单元体属于弯曲下沉带，若0<ε1 t<εc时，该单元体属于裂隙带，若ε1 t≥εc时，该单元体属于垮落带；其中，εc根据实测的地表下沉量与数值模拟结果对比分析确定。
(5)根据垮落带单元体的碎胀系数(BF)判断采空区是否已被完全充填。
所述垮落带单元体的碎胀系数(BF)的表达式为：
其中，hc为煤层厚度，单位为m；hcr为垮落带岩层的厚度，单位为m；采用碎胀系数来确定采空区是否被垮落带岩体充填满，当BF≥1.43时，说明采空区未被充填满；当BF<1.43时，说明采空区已经被充填满。
(6)将垮落带的单元模型由空模型改为双屈服模型，同时弱化裂隙带范围内岩体的材料属性。
材料属性的具体过程包括以下步骤：
A：数值模拟计算中垮落带充填体的相关力学参数的计算主要采用下式：
其中，E0为初始切线模量，MPa；σc为垮落带岩体的单轴抗压强度，MPa；b为垮落带岩体的碎胀系数，即b＝BF。
在FLAC3D软件中，现将垮落带单元体变为空模型，即利用Fish语言中的z_model(_z)＝'null'实现；然后将将空模型行改为双屈服Double-Yield模型，即利用Fish语言中的model doubleyield range group cave_in，并对双屈服模型进行重新赋值。其中的cap压力与单元体塑性体积应变之间的关系通过理论计算可得，但是其余的相关参数还无法确定，诸如密度、体积模量、剪切模量，等等，鉴于无法实验室验证的基础上，可由数值模拟实现——建立一个简单的单元体，单元体的尺寸为1m×1m×1m，进行单轴压缩试验，在单元体上表面施加一固定速度并固定边界。通过不断调整密度、体积模量、剪切模量等物理参量使其与Salamon理论计算曲线吻合。
其中，Salamon表达式为：
其中，σ为当受到侧向约束时加载在松散岩体上的压应力；E0为松散岩体的初始切向应力；εv为体积应变；为最大体积应变；其中，表达式为：
其中，BF为垮落带单元体的碎胀系数。
B：数值模拟计算中裂隙带岩体可等效为弹性各向同性材料，主要改变的是弹性模量和泊松比的降低。研究结果表明，一般裂隙带岩体的弹性模量为500MPa，泊松比为0.3。
(7)重复执行上述过程，直至工作面开采结束。
针对某矿的3108工作面采掘过程中强矿震显现频发情况进行分析说明。
一种基于采空区压实效应的工作面开采全过程数值模拟方法，包括以下步骤：
(1)以某矿3108工作面为模型基础，在FLAC3D软件中建立三维数值计算模型，模型的长、宽、高尺寸为500m×500m×396m。根据现场地应力实测结果，模型中垂直应力和水平应力分别为σzz＝2MPa，σyy＝16MPa，σxx＝24MPa。收集到如表1中的工作面相关地质信息，模型采用应***化模型进行计算。
表1
(2)以某矿3108工作面为模型基础，建立三维数值计算模型，模型的长、宽、高尺寸为500m×500m×396m，如图2所示。根据现场地应力实测结果，模型中垂直应力和水平应力分别为σzz＝2MPa，σyy＝16MPa，σxx＝24MPa，模型采用应***化模型进行计算。
(3)根据建立的模型大小和数量，确定煤层的最小网格大小为5m，因此开挖步距定为一个网格的大小，即5m；当工作面开挖前6个网格时(即开挖30m)，模型运算平衡且最大垂直位移均小于8m(煤层厚度)，当开挖第7个网格(即开挖35m)时，模型运算至最大垂直位移均达到8m时停止运算，因此，可推断该模型的顶板初次来压步距范围为30～35m，下面的计算取初次来压步距为35m。使用上述方法同样进行顶板周期来压步距的运算，以35米处为起点，继续开挖第8个网格(即开挖5m)，模型运算平衡且最大垂直位移均小于8m(煤层厚度)，当开挖第11个网格(即开挖20m)时，模型运算至最大垂直位移均达到8m时停止运算，根据计算结果推测周期来压步距范围为15～20m，下面的计算取周期来压步距20m，依次开挖20m并计算垮落带和裂隙带的范围。
(4)根据实测的地表下沉量与数值模拟结果对比分析确定顶板临界垮落塑性应变值εc＝0.05，因此确定当顶板某一单元体塑性应变值为ε1 t＝0时，该单元体属于弯曲下沉带；当顶板某一单元体塑性应变值为0<ε1 t<0.05时，该单元体属于裂隙带；顶板某一单元体塑性应变值为ε1 t≥0.05时，该单元体属于垮落带。
(5)垮落带单元体的碎胀系数(BF)具体如下：
其中，hc为煤层厚度，单位为m；hcr为垮落带岩层的厚度，单位为m。
根据数值模拟结果可知，hc＝8m，hcr＝30.98m，则BF＝1.26。
(6)数值模拟计算中垮落带充填体的相关力学参数的计算主要采用下式：
其中，E0为初始切线模量，MPa；σc为垮落带岩体的单轴抗压强度，MPa；b为垮落带岩体的碎胀系数，即b＝BF。
垮落带参数拾取的验证：双屈服Double-Yield模型中的cap压力与单元体塑性体积应变之间的关系通过理论计算可得，但是其余的相关参数还无法确定，诸如密度、体积模量、剪切模量，等等，鉴于无法实验室验证的基础上，可由数值模拟实现。建立一个简单的单元体，单元体的尺寸为1m×1m×1m，进行单轴压缩试验，在单元体上表面施加一固定速度并固定边界。通过不断调整密度、体积模量、剪切模量等物理参量使其与Salamon理论计算曲线吻合，如图3所示。最终得到符合数值模拟要求垮落带岩体的的力学参数如表2和表3所示。
表2
表3
裂隙带岩体的材料可等效为弹性各向同性材料，弹性模型参数主要是根据GSI分类进行计算，本发明实施例将裂隙带参数也进行相应的改变，主要改变的是弹性模量和泊松比的降低，本文中裂隙带岩层的弹性模量为160MPa、泊松比为0.3，则相应的体积模量和剪切模量分别为133.33MPa、61.54Mpa，即FLAC3D软件中的Fish语言：z_prop(_z,"bulk")和z_prop(_z,"shear")，具体的可根据矿井工作面垮落带和裂隙带实测的物理力学属性进行修正。
(7)重复执行步骤(1)～(6)，直至工作面开采结束，其工作面开采过程中覆岩垮落带和裂隙带演化特征如图4和图5所示。当工作面首先开挖5m后，覆岩未产生垮落带，但顶板开始破裂，形成裂隙带，高度为8.86m；当工作面继续继续开挖10m后，顶板岩层开始垮落下沉，垮落带高度为2.87m，裂隙带继续往上发展，且裂隙带高度增加至15.35m；当工作面继续进行开挖进程后，覆岩垮落带和裂隙带高度增大，直至达到顶板初次来压(工作面开挖35m)，在此过程中垮落带高度依次为8.86m、13.80m、18.23m、26.73m、30.98m；裂隙带高度依次为18.10m、27.96m、30.32m、34.01m、44.21m。当工作面继续开挖后，可以发现尽管开采后期存在一些奇异点，但考虑到破断、垮落岩层的连贯性，可认为垮落带高度基本保持在30.98m不变；而裂隙带依旧往上发展，直至工作面开挖至95m左右，裂隙带高度保持恒定不变，此时处于工作面见方阶段；再往后工作面上方的垮落带和裂隙带皆已发育至最大高度，垂直方向上无变化，仅仅只是在走向方向上随工作面的推进范围的进一步扩大而已。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。
Claims (8)
1.基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
包括如下步骤：
(1)收集待评估的工作面的相关地质信息；
(2)根据所述地质信息建立数值计算模型；
(3)根据开挖步距和顶板下沉量确定工作面初次来压步距和周期来压步距；
(4)根据计算过程中单元体塑性应变的变化确定工作面上方垮落带和裂隙带的范围，并将垮落带单元体定位为空模型；
(5)根据垮落带单元体的碎胀系数判断采空区是否已被完全充填；
(6)将垮落带的单元模型由空模型改为双屈服模型，同时弱化裂隙带范围内岩体的材料属性；
(7)重复执行上述过程，直至工作面开采结束。
2.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(1)中，相关地质信息包括：岩层的岩性、煤岩体的物理力学属性、工作面埋深、工作面长度以及原岩应力场；
煤岩体的物理力学属性包括煤岩体厚度、弹性模量、泊松比、密度、内摩擦角、内聚力和单轴抗拉强度。
3.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(2)中，建立数值计算模型的方法如下：
建立模型网格，模块参数设置，边界条件与初始条件设置，模型初始平衡。
4.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(3)中，确定工作面初次来压和周期来压步距的方法包括如下步骤：
(31)根据数值模型网格的大小确定开挖步距，设定开挖步距为一个网格；
(32)每运行100步，进行最大不平衡力收敛确定，若最大平衡比率小于1e-5，说明模型受力平衡，可进行进一步的开挖计算；若最大平衡比率大于1e-5，说明模型受力不平衡，此时进行下一步(33)的顶板下沉量判断；
(33)每进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移大于等于开挖煤层的厚度，则说明垮落充分，此时选取该开挖步距内任一点的开挖值为初次来压步距，若最大垂直位移小于开挖煤层的厚度，则继续开挖；
(34)以垮落充分的开挖步距为起点，每继续进行一个开挖步距，根据顶板下沉量，判断顶板是否垮落充分；若最大垂直位移大于等于开挖煤层的厚度，则选取当前开挖步距内任一点的开挖值为周期来压步距，若最大垂直位移小于开挖煤层的厚度，则继续开挖。
5.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(4)中，确定工作面上方垮落带和裂隙带的范围的方法包括如下方法：
基于塑性理论，最大主应变ε1 t表达式如下：
ε1 t＝ε1 e+ε1 p；
其中，ε1 e为最大主弹性应变，ε1 p为最大主塑性应变，εc为顶板岩层单元体垮落前的最大拉应变；若单元体的塑性应变值ε1 t＝0时，该单元体属于弯曲下沉带，若0<ε1 t<εc时，该单元体属于裂隙带，若ε1 t≥εc时，该单元体属于垮落带；其中，εc根据实测的地表下沉量与数值模拟结果对比分析确定。
6.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(5)中，垮落带单元体的碎胀系数的表达式为：
其中，hc为煤层厚度，单位为m；hcr为垮落带岩层的厚度，单位为m；
采用碎胀系数来确定采空区是否被垮落带岩体充填满，当BF≥1.43时，说明采空区未被充填满；当BF<1.43时，说明采空区已经被充填满。
7.根据权利要求1所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(6)中，材料属性的计算方法包括如下内容：
数值模拟计算中垮落带充填体的相关力学参数的计算主要采用下式：
其中，E0为初始切线模量，MPa；σc为垮落带岩体的单轴抗压强度，MPa；b为垮落带岩体的碎胀系数；
建立一个尺寸为1m×1m×1m的单元体进行单轴压缩试验，在单元体上表面施加一固定速度并固定边界；通过不断调整密度、体积模量、剪切模量等物理参量使其与Salamon理论计算曲线吻合，最终得到符合数值模拟要求垮落带岩体的的力学参数；
Salamon表达式为：
其中，σ为当受到侧向约束时加载在松散岩体上的压应力；E0为松散岩体的初始切向应力；εv为体积应变；为最大体积应变；其中，表达式为：
其中，BF为垮落带单元体的碎胀系数；
通过FLAC3D软件将空模型行改为双屈服模型，然后利用上述步骤计算所得的相应力学参数对双屈服模型进行重新赋值；
将裂隙带岩体可等效为弹性各向同性材料，然后降低其弹性模量和泊松比。
8.根据权利要求7所述的基于采空区压实效应的工作面开采全过程数值模拟方法，其特征在于：
步骤(6)中，一般裂隙带岩体的弹性模量为500MPa，泊松比为0.3，相应的单元体体积模量和剪切模量分别为133.36MPa和61.53MPa。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810505702.3A CN108763725A (zh) | 2018-05-24 | 2018-05-24 | 基于采空区压实效应的工作面开采全过程数值模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810505702.3A CN108763725A (zh) | 2018-05-24 | 2018-05-24 | 基于采空区压实效应的工作面开采全过程数值模拟方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108763725A true CN108763725A (zh) | 2018-11-06 |
Family
ID=64005286
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810505702.3A Pending CN108763725A (zh) | 2018-05-24 | 2018-05-24 | 基于采空区压实效应的工作面开采全过程数值模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108763725A (zh) |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109214135A (zh) * | 2018-11-08 | 2019-01-15 | 淮阴工学院 | 一种充填协同垮落式综采面过渡区域局部强矿压控制方法 |
| CN109583116A (zh) * | 2018-12-10 | 2019-04-05 | 西南交通大学 | 基于多源实测信息的地下工程初始应力场动态反演方法 |
| CN109630088A (zh) * | 2018-12-05 | 2019-04-16 | 中煤科工集团重庆研究院有限公司 | 一种高位钻孔布置位置的确定方法 |
| CN109681272A (zh) * | 2018-12-24 | 2019-04-26 | 江西理工大学 | 胶结矿柱支护下金属矿采空区覆岩失稳突变判别方法 |
| CN109783947A (zh) * | 2019-01-21 | 2019-05-21 | 南京大学 | 一种采水型地裂缝数值模拟方法 |
| CN109977453A (zh) * | 2019-01-15 | 2019-07-05 | 河北工程大学 | 固体充填液压支架工作阻力设计方法 |
| CN110362945A (zh) * | 2019-07-22 | 2019-10-22 | 安徽理工大学 | 一种基于flac3d内置fish语言的采动塑性区体积确定方法 |
| CN111428357A (zh) * | 2020-03-20 | 2020-07-17 | 山西工程技术学院 | 基于覆岩剩余自由空间高度的地表最大下沉值确定方法 |
| CN115753442A (zh) * | 2022-11-01 | 2023-03-07 | 中国地质大学（北京） | 适用于煤矿采空区覆岩及地表变形的数值模拟方法及装置 |
| CN115840901A (zh) * | 2022-11-03 | 2023-03-24 | 中国矿业大学 | 一种基于支架阻力连续性分类的工作面来压特征分析方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107942380A (zh) * | 2017-10-23 | 2018-04-20 | 太原理工大学 | 一种考虑垮落带采空区的计算机数值模拟方法 |
-
2018
- 2018-05-24 CN CN201810505702.3A patent/CN108763725A/zh active Pending
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107942380A (zh) * | 2017-10-23 | 2018-04-20 | 太原理工大学 | 一种考虑垮落带采空区的计算机数值模拟方法 |
Non-Patent Citations (2)
| Title |
|---|
| MAHDISHABANIMASHCOOL等: "Numerical modelling of longwall mining and stability analysis of the gates in a coal mine", 《INTERNATIONAL JOURNAL OF ROCK MECHANICS & MINING SCIENCES》 * |
| 朱广安: "深地超应力作用效应及孤岛工作面整体冲击失稳机理研究", 《中国优秀博硕士学位论文全文数据库(博士)工程科技Ⅰ辑》 * |
Cited By (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109214135A (zh) * | 2018-11-08 | 2019-01-15 | 淮阴工学院 | 一种充填协同垮落式综采面过渡区域局部强矿压控制方法 |
| CN109214135B (zh) * | 2018-11-08 | 2022-07-15 | 淮阴工学院 | 一种充填协同垮落式综采面过渡区域局部强矿压控制方法 |
| CN109630088A (zh) * | 2018-12-05 | 2019-04-16 | 中煤科工集团重庆研究院有限公司 | 一种高位钻孔布置位置的确定方法 |
| CN109583116A (zh) * | 2018-12-10 | 2019-04-05 | 西南交通大学 | 基于多源实测信息的地下工程初始应力场动态反演方法 |
| CN109681272A (zh) * | 2018-12-24 | 2019-04-26 | 江西理工大学 | 胶结矿柱支护下金属矿采空区覆岩失稳突变判别方法 |
| CN109977453B (zh) * | 2019-01-15 | 2023-04-18 | 河北工程大学 | 固体充填液压支架工作阻力设计方法 |
| CN109977453A (zh) * | 2019-01-15 | 2019-07-05 | 河北工程大学 | 固体充填液压支架工作阻力设计方法 |
| CN109783947B (zh) * | 2019-01-21 | 2020-12-25 | 南京大学 | 一种采水型地裂缝数值模拟方法 |
| CN109783947A (zh) * | 2019-01-21 | 2019-05-21 | 南京大学 | 一种采水型地裂缝数值模拟方法 |
| CN110362945A (zh) * | 2019-07-22 | 2019-10-22 | 安徽理工大学 | 一种基于flac3d内置fish语言的采动塑性区体积确定方法 |
| CN111428357A (zh) * | 2020-03-20 | 2020-07-17 | 山西工程技术学院 | 基于覆岩剩余自由空间高度的地表最大下沉值确定方法 |
| CN111428357B (zh) * | 2020-03-20 | 2023-03-28 | 山西工程技术学院 | 基于覆岩剩余自由空间高度的地表最大下沉值确定方法 |
| CN115753442A (zh) * | 2022-11-01 | 2023-03-07 | 中国地质大学（北京） | 适用于煤矿采空区覆岩及地表变形的数值模拟方法及装置 |
| CN115753442B (zh) * | 2022-11-01 | 2023-09-05 | 中国地质大学（北京） | 适用于煤矿采空区覆岩及地表变形的数值模拟方法及装置 |
| CN115840901A (zh) * | 2022-11-03 | 2023-03-24 | 中国矿业大学 | 一种基于支架阻力连续性分类的工作面来压特征分析方法 |
| CN115840901B (zh) * | 2022-11-03 | 2023-09-05 | 中国矿业大学 | 一种基于支架阻力连续性分类的工作面来压特征分析方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108763725A (zh) | 基于采空区压实效应的工作面开采全过程数值模拟方法 | |
| Li et al. | Large scale three-dimensional seepage analysis model test and numerical simulation research on undersea tunnel | |
| CN104458309B (zh) | 一种用于物理模拟实验中的相似材料配比确定方法 | |
| CN108868770B (zh) | 一种充填开采岩层位态精准控制设计方法 | |
| CN102866241B (zh) | 三向加载大型三维相似模拟试验方法 | |
| CN103983742B (zh) | 煤层覆岩破断煤岩体瓦斯运移及抽采实验*** | |
| CN103344491B (zh) | 一种基于静载和***荷载共同作用下巷道冲击地压的模拟方法 | |
| CN105550441B (zh) | 一种基于连续介质的工程岩体破裂劣化数值模拟方法 | |
| Lee et al. | Rock engineering in underground energy storage in Korea | |
| Detournay et al. | FLAC and numerical modeling in geomechanics | |
| Al Heib et al. | On the use of 1g physical models for ground movements and soil-structure interaction problems | |
| CN102879284B (zh) | 三向加载大型三维相似模拟试验试件箱 | |
| CN109492262A (zh) | 一种利用数值模拟分析非均匀分布裂隙巷道稳定性的方法 | |
| CN111595703A (zh) | 一种基于模型试验的节理岩质边坡***失稳规律的测试方法 | |
| CN113803066A (zh) | 无煤柱自成巷平衡开采设计方法及*** | |
| CN102879550B (zh) | 三向加载大型三维相似模拟试验载荷模拟方法 | |
| CN115903078A (zh) | 水力压裂层位及其确定方法、水力压裂位置的确定方法 | |
| Sun et al. | Physical model experiment and numerical analysis on innovative gob-side entry retaining with thick and hard roofs | |
| Arasteh et al. | Discontinuous modeling of roof strata caving in a mechanized longwall mine in tabas coal mine | |
| CN111695790A (zh) | 一种保安矿柱开采方法 | |
| CN205353043U (zh) | 一种充填注浆试验装置 | |
| Wang et al. | Shaking table tests and numerical analysis on the seismic response of karst-crossing socketed piles in dry sandy soil foundation | |
| Moussaei et al. | Physical modeling of soil arching around shallow tunnels in sandy grounds | |
| Kalenchuk et al. | Three-dimensional numerical simulations of the Downie Slide to test the influence of shear surface geometry and heterogeneous shear zone stiffness | |
| CN108664743A (zh) | 一种基于覆岩演化特征的煤矿动载显现数值模拟方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181106 |