D010.13374f.isn0W3x.20l.02.001 第33卷第2期 北京科技大学学报 Vo133N92 2011年2月 Journal ofUniversity of Science and Technobgy Bejjing Feb 2011 采空区坚硬顶板流变破断力学分析 王金安” 李大钟”尚新春D 1)北京科技大学土木与环境工程学院。北京1000832)香港理工大学土木学院。香港 通信作者，Em时w@ub业m 摘要考虑岩体的流变特性，运用弹性一黏弹性对应准则及【即数值逆变换对矿柱支撑的采空区顶板岩层变形进行了 分析，建立了采空区顶板挠度随时间的变化关系，揭示了采空区项板位移随时间变化的特性.研究表明，在考虑岩体流变性质 的情况下，采空区顶板下沉位移随时间延长而增大，在给定比较长的时间跨度，即使坚硬岩体项板位移仍会达到相当大的量 值.进而引起采空区项板破断.通过对邢台石膏矿塌陷案例分析，验证了计算方法的可行性。 关键词矿业工程：顶板控制：岩石力学；蠕变：数学模型 分类号TD322 M echanics anaysis on creep fracture of strong roof strata above m ned-out area WANG Jinan)LIDa_ong 2)SHANG Xin chuni) 1)ShoolofCivil and Envirormenl Engineering University of Science and Technokgy Beijing Beijng 100083 Chna 2)ShoolofCivil Engneerng Hong Kang Pol technic University HongK ang China Correspand ng author Email wp@ustb edu cn ABSTRACT The depmation of roof strata suppored by pillar n mined out area was analyzed n consderation of he theological Prperties of rock masses and by utilizing he elstic viscoe lastic correspondence principle and mumerical laplace nverson The roof defecton as a fnction of tme was established based on which the tme dependentbehavpr of the roofd isplacementwas elaborated It is shown that when rock creep properties are taken np account pr both he pillars and roof the roofd isplcement ncreases with pro longing tme as expected Therepre given a sufficient perod of tie the roof depmaton can reach a consderab le magnitude even pr had rock masses and even ualy lead p the collapse of roof stat As a case the estination of filre tine was camried out prX ngtai Gypsum Mne in Ch ina in wh ich the feasibility of he poposedmethod was verifed by acceptable agreement beween the estiated da ta and actual field observatons KEY WORDS minng engineering mine roof contrl rock mechan ics creep mathematical models 采空区顶板的大规模塌陷将造成巨大的财产损 某一输出量进行预测，采用以案例统计为基础的数 失和人员伤亡).采用合理的方法估计采空区稳定 学或人工智能的方法一直以来被认为有广阔的应用 时间对矿山开采及采空区的处理具有十分重要意 前景，有学者将此类方法应用于采空区稳定性时间 义，由于采空区稳定性特性随时间变化影响因素复 评价中2-.但是，由于我国对采空区稳定影响指标 杂，包括顶板、矿柱流变风化，侵蚀以及应力剥落 统计缺乏以及统计数据可靠度不高，所以目前采用 等，一直以来未形成一套公认的采空区稳定时间预 此类方法缺乏必要的数据基础.通过建立力学模 测体系.目前在此方面研究大致可分为两种方法： 型，王金安等“提出采用类似Wke弹性地基梁 一种是抛开采空区力学特性、破坏机理的数学及人 假设，将矿柱等效为均匀弹簧，顶板抽象为弹性薄板 工智能方法：另一方法则是从风化、侵蚀和岩体流变 从而得出矿柱总面积和顶板力学状态之间的关系， 等方面入手，建立数学和力学模型来解释空区稳定 这实质上是考虑风化和剥落减小矿柱面积对采空区 性时间相依性，从而来预测稳定时间.对复杂系统 稳定性状态的影响，但文中未涉及如何将矿柱面积 收稿日期：2010-04-20 基金项目：国家高技术研究发展计划资助项目(N92008A062104片国家重点基础研究发展计划资助项目(NQ2010CB731500)
第 33卷 第 2期 2011年 2月 北 京 科 技 大 学 学 报 JournalofUniversityofScienceandTechnologyBeijing Vol.33 No.2 Feb.2011 采空区坚硬顶板流变破断力学分析 王金安 1) 李大钟 1, 2) 尚新春 1) 1)北京科技大学土木与环境工程学院, 北京 100083 2)香港理工大学土木学院, 香港 通信作者, E-mail:wja@ustb.edu.cn 摘 要 考虑岩体的流变特性, 运用弹性--黏弹性对应准则及 Laplace数值逆变换对矿柱支撑的采空区顶板岩层变形进行了 分析, 建立了采空区顶板挠度随时间的变化关系, 揭示了采空区顶板位移随时间变化的特性.研究表明, 在考虑岩体流变性质 的情况下, 采空区顶板下沉位移随时间延长而增大, 在给定比较长的时间跨度, 即使坚硬岩体顶板位移仍会达到相当大的量 值, 进而引起采空区顶板破断.通过对邢台石膏矿塌陷案例分析, 验证了计算方法的可行性. 关键词 矿业工程;顶板控制;岩石力学;蠕变;数学模型 分类号 TD322 Mechanicsanalysisoncreepfractureofstrongroofstrataabovemined-outarea WANGJin-an1) , LIDa-zhong1, 2) , SHANGXin-chun1) 1)SchoolofCivilandEnvironmentalEngineering, UniversityofScienceandTechnologyBeijing, Beijing, 100083, China 2)SchoolofCivilEngineering, HongKongPolytechnicUniversity, HongKong, China Correspondingauthor, E-mail:wja@ustb.edu.cn ABSTRACT Thedeformationofroofstratasupportedbypillarsinminedoutareawasanalyzedinconsiderationoftherheological propertiesofrockmassesandbyutilizingtheelastic-viscoelasticcorrespondenceprincipleandnumericalLaplaceinversion.Theroof deflectionasafunctionoftimewasestablished, basedonwhichthetime-dependentbehavioroftheroofdisplacementwaselaborated. Itisshownthatwhenrockcreeppropertiesaretakenintoaccountforboththepillarsandroof, theroofdisplacementincreaseswithprolongingtimeasexpected.Therefore, givenasufficientperiodoftime, theroofdeformationcanreachaconsiderablemagnitudeevenfor hardrockmassesandeventuallyleadtothecollapseofroofstrata.Asacase, theestimationoffailuretimewascarriedoutforXingtai GypsumMineinChina, inwhichthefeasibilityoftheproposedmethodwasverifiedbyacceptableagreementbetweentheestimateddataandactualfieldobservations. KEYWORDS miningengineering;mineroofcontrol;rockmechanics;creep;mathematicalmodels 收稿日期:2010--04--20 基金项目:国家高技术研究发展计划资助项目(No.2008AA062104);国家重点基础研究发展计划资助项目(No.2010CB731500) 采空区顶板的大规模塌陷将造成巨大的财产损 失和人员伤亡 [ 1] .采用合理的方法估计采空区稳定 时间对矿山开采及采空区的处理具有十分重要意 义 .由于采空区稳定性特性随时间变化影响因素复 杂 ,包括顶板 、矿柱流变, 风化 , 侵蚀以及应力剥落 等 ,一直以来未形成一套公认的采空区稳定时间预 测体系 .目前在此方面研究大致可分为两种方法 : 一种是抛开采空区力学特性 、破坏机理的数学及人 工智能方法 ;另一方法则是从风化、侵蚀和岩体流变 等方面入手 ,建立数学和力学模型来解释空区稳定 性时间相依性, 从而来预测稳定时间.对复杂系统 某一输出量进行预测 ,采用以案例统计为基础的数 学或人工智能的方法一直以来被认为有广阔的应用 前景 ,有学者将此类方法应用于采空区稳定性时间 评价中 [ 2--3] .但是, 由于我国对采空区稳定影响指标 统计缺乏以及统计数据可靠度不高 , 所以目前采用 此类方法缺乏必要的数据基础 .通过建立力学模 型, 王金安等 [ 1, 4] 提出采用类似 Winkler弹性地基梁 假设 ,将矿柱等效为均匀弹簧 ,顶板抽象为弹性薄板 从而得出矿柱总面积和顶板力学状态之间的关系, 这实质上是考虑风化和剥落减小矿柱面积对采空区 稳定性状态的影响 ,但文中未涉及如何将矿柱面积 DOI :10 .13374 /j .issn1001 -053x .2011 .02 .001
第2期 王金安等：采空区坚硬顶板流变破断力学分析 143 减小与采空区稳定时间等关联的讨论.Bwa对 弹性问题的解可以通过对弹性解作Laplacei逆变换 不同时间的矿柱进行了应力测量，回归了风化厚度 得到.该方法严密的论证以及局限性可参见文 同矿柱时间的关系.Auvray等研究了相对湿度对 献[13]. 矿柱强度的衰减.Caste lanza等根据水影响下石 总之，对于静态黏弹性问题，采用弹性黏弹性 膏矿柱单轴压缩试验数据建立风化模型进而对矿柱 对应准则，求解过程可分为：①求解对应的弹性问 稳定时间做出预测. 题，②对弹性解进行Laplace变换；③用松弛模量替 岩土工程中考虑岩土材料的流变特性，对设计 代弹性模量：④进行LPe逆变换便可求得问题 施工均具有鲜明而重要的意义【89，在地下隧洞位 的解. 移时间特性分析中，得到了广泛的应用心四.本文 本文将矿柱受力视为单轴压缩，顶板偏应力应 针对矿柱支撑下坚硬顶板的大面采空区，建立考虑 变关系采用Bu您模型，顶板球应力应变张量关系 流变作用采空区矿柱顶板变形分析的力学模型，揭 采用弹性模型.根据以上分析，先对模型求弹性解. 示变形和顶板破坏随时间的变化特征 1.1模型弹性解 针对顶板较为坚硬而矿柱相对软弱的大面积采 1采空区顶板的黏弹性力学模型 空区，将顶板简化为弹性矩形薄板，矿柱等效成均匀 在应力应变分析中，共需要平衡方程、相容方 分布黏弹性地基，建立如图1模型.其中弹性矩形 程、本构方程及边界条件.对于弹性和黏弹性分析， 平板长度为2宽度为2b(起)和厚度为h顶板 前两类方程并无区别.区分弹性与黏弹性最主要是 岩体的弹性模量为E泊松比为体密度为P和抗 本构关系，通过Laplace变换，在Lpce域上的黏 拉强度σī上层岩土介质对顶板表面的压力为均 弹性问题与弹性问题求解非常类似，通常情况下黏 布载荷. 地表 程盖层 坚使顶板 矿柱 图1采空区示意图()、顶板岩体简化为四边固支的弹性矩形平板()及模型沿一剂面图() Fig I Schem atic illustra tion ofm ined out aea a)smplified elstic rectangular plate for roof stmta (b)and cross sectional pofile on x-zplane (c) 矿柱平均面积为A高度为H假设矿柱是等距 9.顶板在破坏前为固支边界条件： 分布的，其总数目为P弹性假设下采空区顶板的控 w±a=0w±b=0 制方程为 w w ±=0别-=0 (2) D△4w+kw=q (1) 1.2顶板边缘破裂前的位移方程 式中：w为顶板挠度；D为板抗弯刚度，D=ER/ 根据顶板挠度和边缘约束条件，取试函数为 [12(1一2)]，E为顶板弹性模量；k=EAy (4aH,H为矿柱高度；E为随矿柱时间变化的弹 爬(-(-B时 (3) 性模量：为顶板自重与承受外部荷载之和，4-P叶 将其代入式(1)的伽辽金弱形式方程，可得顶板中
第 2期 王金安等:采空区坚硬顶板流变破断力学分析 减小与采空区稳定时间等关联的讨论.Biswas [ 5] 对 不同时间的矿柱进行了应力测量, 回归了风化厚度 同矿柱时间的关系.Auvray等 [ 6]研究了相对湿度对 矿柱强度的衰减 .Castellanza等 [ 7] 根据水影响下石 膏矿柱单轴压缩试验数据建立风化模型进而对矿柱 稳定时间做出预测. 岩土工程中考虑岩土材料的流变特性, 对设计 施工均具有鲜明而重要的意义 [ 8--9] , 在地下隧洞位 移时间特性分析中 ,得到了广泛的应用 [ 10--12] .本文 针对矿柱支撑下坚硬顶板的大面采空区, 建立考虑 流变作用采空区矿柱顶板变形分析的力学模型 ,揭 示变形和顶板破坏随时间的变化特征 . 1 采空区顶板的黏弹性力学模型 在应力应变分析中, 共需要平衡方程 、相容方 程 、本构方程及边界条件.对于弹性和黏弹性分析 , 前两类方程并无区别 .区分弹性与黏弹性最主要是 本构关系 , 通过 Laplace变换, 在 Laplace域上的黏 弹性问题与弹性问题求解非常类似, 通常情况下黏 弹性问题的解可以通过对弹性解作 Laplace逆变换 得到 .该方法严密的论证以及局限性可参见文 献 [ 13] . 总之 ,对于静态黏弹性问题 ,采用弹性 --黏弹性 对应准则 ,求解过程可分为 :①求解对应的弹性问 题;②对弹性解进行 Laplace变换;③用松弛模量替 代弹性模量;④进行 Laplace逆变换便可求得问题 的解 . 本文将矿柱受力视为单轴压缩, 顶板偏应力应 变关系采用 Burgers模型 ,顶板球应力应变张量关系 采用弹性模型 .根据以上分析, 先对模型求弹性解. 1.1 模型弹性解 针对顶板较为坚硬而矿柱相对软弱的大面积采 空区 ,将顶板简化为弹性矩形薄板 ,矿柱等效成均匀 分布黏弹性地基 ,建立如图 1模型 .其中弹性矩形 平板长度为 2a、宽度为 2b(b≤a)和厚度为 h, 顶板 岩体的弹性模量为 Er、泊松比为 ν、体密度为 ρ和抗 拉强度 σT.上层岩土介质对顶板表面的压力为均 布载荷. 图 1 采空区示意图(a)、顶板岩体简化为四边固支的弹性矩形平板(b)及模型沿 x-z剖面图(c) Fig.1 Schematicillustrationofminedoutarea(a), simplifiedelasticrectangularplateforroofstrata(b)andcrosssectionalprofileonx-zplane (c) 矿柱平均面积为 A,高度为 H.假设矿柱是等距 分布的 ,其总数目为 n, 弹性假设下采空区顶板的控 制方程为 [ 4] DΔ 4w+kw=q (1) 式中:w为顶板挠度;D为板抗弯刚度, D=Erh 3 / [ 12(1 -ν 2 )] , Er为 顶板弹性模 量;k =nEpA/ (4abH), H为矿柱高度;Ep为随矿柱时间变化的弹 性模量;q为顶板自重与承受外部荷载之和, q=ρgh+ q0 .顶板在破坏前为固支边界条件 : w x=±a =0, w y=±b =0, w x x=±a =0, w y y=±b =0 (2) 1.2 顶板边缘破裂前的位移方程 根据顶板挠度和边缘约束条件,取试函数为 w= w0 a 4 b 4 (x 2 -a 2 ) 2(y 2 -b 2) 2 (3) 将其代入式 (1)的伽辽金弱形式方程 , 可得顶板中 · 143·
。144 北京科技大学学报 第33卷 面的最大挠度（中心挠度）为 w= 69 (4) E fE+ -T一 式中， 子++ 图2Buer体物理本构模型 Fg2 Bugers physical constitutive molel 1.3顶板边缘破坏之后位移方程 四边铰支的边界条件为 wl士=0wsb=0 寻w a w =0丽±=0 0汉无士n (5) 取试函数为 ww cot xcot y (6) 2a 2b 同理，由迦辽金法求得 ca 图3Bu呕er体蠕变曲线 w- E+FE (7) Fig 3 Creep curve for the Burgersmodel 1~ 式中，- g=π2i1+12 式种，用=餐++号母=装9=9= 192aTB· 1112 2微分形式的线性黏弹性本构关系 根据式(11)和式（12，对于Bus体的 采用线性黏弹性模型，简单应力状态下如单 Laplace变换等效松弛模量为 轴压缩或纯剪)微分形式的材料本构关系可以表 0 Q9=9=9斗93 示为 EP 1+月+B3 (13) Po=Q (8) 21复杂应力状态下微分形式的黏弹性本构模型 式中，P和Q为对时间的线性微分算子，形式如下： 复杂应力状态下，由于黏弹性材料对剪应力 会易Q-宫月 (9) (偏应力)和静水应力响应不同，习惯上将黏弹性本 构关系表示为偏张量和球张量的形式： 对式(8)进行Laplace变换（初值条件为0） B(9=Qd(9 得到 (14) P(今G(习=(B+胃斗93+…十BS)。(9= Boi(=QEi( 式中，、B、Q和Q为对时间的微分算子，同 Q可et9=(9+9斗g3+.+9$)et9 式(8，和份别为应力和应变张量，σ和e分 (10) 别为应力和应变球张量. 式中，为L即ce变量. 对于弹性各向同性材料，本构关系可以表示为： 由式(10)Laplac变换后松弛模量可以表示为 =2Gdj 。Q9=9 (15) (11) eP(可 0i=3Ki 通过式(14入式(15河以得出如下关系： 由Kevm体和Mawe体串联组成的Bugers 体，具有四个可调参数，可以描述第三期蠕变前的蠕 G 2P K19 3B (16) 变曲线，灵活实用，其物理模型及蠕变曲线如 式中，G和K分别为剪切模量和体积模量.弹性模 图2和图3所示. 量Ev与GK有如下关系： Bugers体本构模型见下式2-， E9KG 张器 (17) G+Po+Ro=9e+9e (12)
北 京 科 技 大 学 学 报 第 33卷 面的最大挠度(中心挠度)为 w0 = c0 q c1 Ep +c2 Er 1 -ν 2 (4) 式中, c0 = 441 128 , c1 = nA 2abH , c2 = 3h 3 4 7 a 4 + 4 a 2 b 2 + 7 b 4 . 1.3 顶板边缘破坏之后位移方程 四边铰支的边界条件为 w x=±a =0, w y=±b =0, 2w x 2 x=±a =0, 2 w y 2 y=±b =0 (5) 取试函数为 w=w0cos πx 2a cos πy 2b (6) 同理, 由迦辽金法求得 w0 = c0′q c′1Ep + c2′Er 1 -ν 2 (7) 式中, c0′= 16 π 2 , c′1 = nA 4abH , c2′= π 2h 3 192 1 a 2 + 1 b 2 2 . 2 微分形式的线性黏弹性本构关系 采用线性黏弹性模型, 简单应力状态下 (如单 轴压缩或纯剪 )微分形式的材料本构关系可以表 示为 Pσ=Qε (8) 式中, P和 Q为对时间的线性微分算子 ,形式如下 : P=∑ a r=0 pr r t r, Q=∑ b r=0 qr r t r (9) 对式 (8)进行 Laplace变换 (初值条件为 0) 得到 P (s)σ (s)=(p0 +p1 s+p2 s 2 +… +pas a)σ (s)= Q (s)ε (s)=(q0 +q1s+q2s 2 +… +qbs b)ε (s) (10) 式中, s为 Laplace变量 . 由式(10)Laplace变换后松弛模量可以表示为 σ ε = Q (s) P (s) =sE (s) (11) 由 Kelvin体和 Maxwell体串联组成的 Burgers 体 ,具有四个可调参数 ,可以描述第三期蠕变前的蠕 变曲线 , 灵活实用 [ 14] , 其物理模型及蠕变曲线如 图 2和图 3所示 . Burgers体本构模型见下式 [ 12--17] : σ+p1 σ · +p2σ ·· =q1 ε · +q2 ε ·· (12) 图 2 Burgers体物理本构模型 Fig.2 Burgersphysicalconstitutivemodel 图 3 Burgers体蠕变曲线 Fig.3 CreepcurvefortheBurgersmodel 式中 , p1 = η1 k1 + η1 k2 + η2 k2 , p2 = η1 η2 k1 k2 , q1 =η1 , q2 = η1 η2 k2 . 根据 式 (11)和式 (12), 对于 Burgers体 的 Laplace变换等效松弛模量为 σ ε = Q (s) P (s) =sE (s)= q1s+q2s 2 1 +p1s+p2s 2 (13) 2.1 复杂应力状态下微分形式的黏弹性本构模型 复杂应力状态下, 由于黏弹性材料对剪应力 (偏应力)和静水应力响应不同 ,习惯上将黏弹性本 构关系表示为偏张量和球张量的形式: P1sij(t)=Q1dij(t) P2σii(t)=Q2 εii(t) (14) 式中 , P1 、 P2 、 Q1 和 Q2 为对时间的微分算子 , 同 式(8).sij和 dij分别为应力和应变张量, σii和 εii分 别为应力和应变球张量 . 对于弹性各向同性材料 ,本构关系可以表示为: sij=2Gdij σii =3Kεii (15) 通过式(14)、式 (15)可以得出如下关系: G= 1 2 Q1 P1 , K= 1 3 Q2 P2 (16) 式中 , G和 K分别为剪切模量和体积模量 .弹性模 量 E、ν与 G、K有如下关系 : E= 9KG 3K+G , ν= 3K-2G 6K+2G (17) · 144·
第2期 王金安等：采空区坚硬顶板流变破断力学分析 145 将式(16)中KG与、Q和Q关系代入 69 %(9= E (20) 式(17，得到 FE 3QQ 1- E9-BQ+2RQ 顶板处于复杂应力状态，设顶板材料在静水压 v(9= PQ-BQ (18) 力下为弹性，偏应力下为Bugers黏弹性模型.则 BQ+2PQ 式(18的Ipac变换形式为 P=1+B斗B5Q=9+9,s (9= 3QQ P2=1Q=3K (21) BQ+2RQ 其中，下标表示顶板参数，其他符号物理意义同 PQ-BQ y(9= (19) 前.将式(21)代入式(18)得LPc变换形式，便 BQ+2RQ 可以得到平（和T的表达式. 式中，Q、、Q和B意义同式(14，形式如 矿柱处于简单应力状态（单轴压缩，可根据 式(10，与选取模型有关.例如： 式(13直接求得 (1)Hook体，P(今=1Q(习=3K 9叶9男 (9-1十，叶3 (22) (2)Bugers体，P(=1+P斗？专Q(习= 9斗93 其中，下标表示矿柱参数. 2.2采空区顶板挠度的黏弹性解 将E(.T9和E(替换式(20)中的弹性 对式(4进行La即pac变换因为弹性解同时间 参数，便可得到Laplacet域上的顶板中心的挠度 无关，则有 函数： 69 %(习= (23) %斗9%5SY9+99(9斗45+6K+6K斗6K四8) 1+R,斗R,3(1+B斗B5)(29+295+3K+3K码+3KB3】 原则上讲，对式(23进行Laplace逆变换便可 得到关于时间的顶板位移函数.然而，式(23的 解析LPac逆变换非常困难，甚至不可能.本文对 (25) 式(23求数值Laplace逆变换. Dub1与Cum进而提出了两个计算公式 类似地，将式(23)中的6、和替换为式(7) 中的&、和弹性参数，然后进行Laplace逆变换， 6(-2言时4 便可得到铰支条件下顶板位移随时间的变化关系. (26) 3求解Lap lace数值逆变换 由于数值反演是病态问题，尤其当时间相对较 (27) 大时.目前己提出了许多不同的反演手段，文献 【I8对LPac数值逆总结为四类方法，并对其特 式(25)~(27)中：ST为计算参数，且N为级 点进行了分析，其中Foure级数展开法在黏弹性问 数截取项数；为虚数单位.运用Laplace数值逆求 题的Laplace反演中适应范围较广.下面对Fourjer 解过程中，参数选取对结果有较大的影响，文 级数法求解LaPc逆的方法进行简单叙述. 献[19对Fairjer级数法的参数选取进行了一定的 Lpc逆变换为 讨论.尽管如此，笔者建议对同一问题采用至少两 千射 种不同的数值反演方法来确认解的一致性可参见 (24) NAG数学库，关于Laplace数值逆的详细讨论可 Fou rier级数展开法将求解f〔D的问题转化为 参见文献[13]. 一个广义积分的问题.进而导出原函数的Fourjer 级数表达式.依此Dubner等给出了〔的计算 4案例分析 格式 2005年11月6日邢台尚旺庄石膏矿区发生特
第 2期 王金安等:采空区坚硬顶板流变破断力学分析 将式 (16)中 K、 G与 P1 、 P2 、 Q1 和 Q2 关系代入 式 (17), 得到 E(t)= 3Q1 Q2 P2Q1 +2P1 Q2 , ν(t)= P1 Q2 -P2 Q1 P2Q1 +2P1Q2 (18) 式 (18)的 Laplace变换形式为 sE(s)= 3Q 1 Q 2 P 2Q 1 +2P 1 Q 2 , ν(s)= P 1 Q 2 -P 2 Q 1 P 2Q 1 +2P 1 Q 2 (19) 式中, Q 1 、 P 1 、 Q 2 和 P 2 意 义同式 (14), 形式如 式 (10), 与选取模型有关.例如 : (1)Hook体, P (s)=1, Q (s)=3K; (2)Burgers体 , P (s)=1 +p1s+p2 s 2 , Q (s)= q1s+q2s 2 . 2.2 采空区顶板挠度的黏弹性解 对式(4)进行 Laplace变换, 因为弹性解同时间 无关, 则有 w 0 (s)= c0q sc1 Ep + c2Er 1 -ν 2 (20) 顶板处于复杂应力状态, 设顶板材料在静水压 力下为弹性,偏应力下为 Burgers黏弹性模型.则 Pr1 =1 +pr1 s+pr2 s 2 , Qr1 =qr1 +qr2s, Pr2 =1, Qr2 =3K (21) 其中 ,下标 r表示顶板参数 , 其他符号物理意义同 前.将式(21)代入式 (18)得 Laplace变换形式 , 便 可以得到 sE r(s)和 v (s)的表达式 . 矿柱处于简单应力状态 (单轴压缩 ), 可根据 式(13)直接求得 sE p(s)= qp1s+qp2s 2 1 +pp1s+pp2s 2 (22) 其中 ,下标 p表示矿柱参数 . 将 sE r(s)、v (s)和 sE p(s)替换式 (20)中的弹性 参数 , 便可得到 Laplace域上的顶板中心的挠度 函数 : w 0(s)= c0 q sc1 qp1s+qp2s 2 1 +pp1s+pp2s 2 + c2s(qr1 +qr2 s)(qr1 s+qr2 s 2 +6K+6Kpr1s+6Kpr2 s 2) (1 +pr1s+pr2s 2 )(2qr1s+2qr2s 2 +3K+3Kpr1s+3Kpr2 s 2) (23) 原则上讲 , 对式 (23)进行 Laplace逆变换便可 得到关于时间 t的顶板位移函数 .然而, 式 (23)的 解析 Laplace逆变换非常困难 ,甚至不可能.本文对 式 (23)求数值 Laplace逆变换 . 类似地 ,将式 (23)中的 c0 、c1和 c2替换为式 (7) 中的 c0′、c′1和 c2′弹性参数, 然后进行 Laplace逆变换 , 便可得到铰支条件下顶板位移随时间的变化关系. 3 求解 Laplace数值逆变换 由于数值反演是病态问题, 尤其当时间相对较 大时 .目前已提出了许多不同的反演手段, 文献 [ 18]对 Laplace数值逆总结为四类方法, 并对其特 点进行了分析,其中 Fourier级数展开法在黏弹性问 题的 Laplace反演中适应范围较广.下面对 Fourier 级数法求解 Laplace逆的方法进行简单叙述 . Laplace逆变换为 f(t)= 1 2πi∫ c+iω c-iω F(s)e stds (24) Fourier级数展开法将求解 f(t)的问题转化为 一个广义积分的问题 .进而导出原函数的 Fourier 级数表达式 .依此 Dubner等 [ 15] 给出了 f(t)的计算 格式 f(Ⅰ )(t)= 2e ct T 1 2 F(c)ReF c+ kπ T icos kπt T (25) Durbin [ 16] 与 Crump [ 17]进而提出了两个计算公式 f(Ⅱ )(t)= 2e ct T ∑ N k=1 ImF c+ kπ T isin kπt T (26) f(t)= e ct T 1 2 F(c)+∑ N k=1 ReF c+ kπ T i - ImF c+ kπ T isin kπt T (27) 式(25)～ (27)中:c、T为计算参数 ,且 T≥t;N为级 数截取项数;i为虚数单位 .运用 Laplace数值逆求 解过程 中, 参数 选取对结 果有较 大的影 响, 文 献[ 19]对 Fourier级数法的参数选取进行了一定的 讨论 .尽管如此 ,笔者建议对同一问题采用至少两 种不同的数值反演方法来确认解的一致性 (可参见 NAG数学库 ).关于 Laplace数值逆的详细讨论可 参见文献 [ 13] . 4 案例分析 2005年 11月 6日邢台尚旺庄石膏矿区发生特 · 145·
。146 北京科技大学学报 第33卷 别重大坍塌事故.事故造成37人死亡井下16 尚汪庄石膏矿区矿体顶板埋深138~162甲底 人，地面21人，88间房屋倒塌，8个竖井严重开裂 板埋深150~216四矿体呈透镜状，走向近南北，倾 变形受损.地表塌陷面积5.3×10塌陷区呈 向东，倾角3°~10°：矿体走向长1200四宽300~ 300m210m椭圆形：坍塌体积24.3×10m.地表 700四矿石类型为角砾状石膏矿石，矿层厚3.96~ 移动面积245X10.地面最大倾斜95mm/m 71.7四平均厚27.4四矿体赋存在奥陶系中统下 (约7°)：最大错动量1.5四塌陷区中部最大下沉 马家沟组第二段中下部，如图5所示.矿层直接顶 8.0m(图4). 板为含泥质很高的钙质黏土岩，厚0.2~1.49间接 顶板为角砾状泥灰岩，矿体顶面距第四纪底面距离 矿山边界 为10.8~60.6四矿体直接底板为含黏土质很高的 钙质黏土岩，厚0.5~4.6甲间接底板为薄层状泥质 降区域 白云岩，部分为角砾状泥质白云岩，矿层底板距奥陶 系下统3.6~19.3四该区主要岩石力学参数见 陷区 表1. 柱状 层厚 岩性 80-170m 结土，砂质黏土和可岩 00 10.8- 角砾状泥灰岩 60.6m 02-1.4m 钙质黏土岩 396-711m 274m 石舟矿 05-4.6m 含黏土质很高的钙质结上岩 3.1-14.7m 薄层状泥质白云岩，部分为 角砾状泥质白云岩 图4邢台县尚旺庄石膏矿区矿井分布 图5矿区地层柱状图 Fg 4 Distrbution of mines in the gypam filed at Shangwang vil Fg 5 Geokgical pofile in them ine area lage X nmiCount in north Chima 表1矿区主要矿岩力学参数 T able Mechan ical param eters of rocks in them nes 岩层 弹性模量MPa 泊松比 容重/(kNr3) 内聚力MPa 内摩擦角度（°） 抗拉强度MPa 表土 500 0.30 1&0 08 20.0 0.05 第四纪沉积岩 8000 023 21.0 032 21.5 026 灰岩 500000 025 240 840 400 470 石普 70000 030 23.6 540 38.6 280 河北省邢台县尚旺庄石膏矿区的五个矿先后建 横截面积占采空区面积的比率为=39.16%，矿柱 于1984一1998年期间，有第二石膏矿、太行石膏矿、 的高度=8四根据流变试验结果，运用遗传算法 林旺石膏矿、邢燕石膏矿和康立石膏矿（原恒昌石 编程进行参数识别心刘，得到矿柱和顶板的流变模 膏矿五座矿山.塌陷区发生在康立矿和林旺矿之 型参数分别为7n=6.41×10MPah72=4.50X 间，塌陷区周边是各矿井的安全矿柱（图4）.将塌 10°MPhk=27.7×10MPk=2.95× 陷区所对应的采空区近似成矩形区域则坚硬顶板 10MP?顶板流变参数为7，=1.03×10MPah 长度为2280四宽度为2b=180四矿柱群的总 7。=468×102 MPa h k=50X103MP9k=
北 京 科 技 大 学 学 报 第 33卷 别重大坍塌事故 [ 1] .事故造成 37 人死亡 (井下 16 人 、地面 21人), 88间房屋倒塌, 8个竖井严重开裂 变形受损.地表塌陷面积 5.3 ×10 4 m 2 , 塌陷区呈 300m×210m椭圆形 ;坍塌体积 24.3 ×10 4 m 3 .地表 移动面积 24.5 ×10 4 m 2.地面最大倾斜 95 mm/m (约 7°);最大错动量 1.5 m;塌陷区中部最大下沉 8.0m(图 4). 图 4 邢台县尚旺庄石膏矿区矿井分布 Fig.4 DistributionofminesinthegypsumfiledatShangwangvillage, XintaiCountyinnorthChina 尚汪庄石膏矿区矿体顶板埋深 138 ～ 162 m, 底 板埋深 150 ～ 216 m;矿体呈透镜状 ,走向近南北, 倾 向东 ,倾角 3°～ 10°;矿体走向长 1 200 m, 宽 300 ～ 700 m,矿石类型为角砾状石膏矿石, 矿层厚 3.96 ～ 71.7 m,平均厚 27.4 m.矿体赋存在奥陶系中统下 马家沟组第二段中下部, 如图 5所示 .矿层直接顶 板为含泥质很高的钙质黏土岩,厚 0.2 ～ 1.4 m;间接 顶板为角砾状泥灰岩 ,矿体顶面距第四纪底面距离 为 10.8 ～ 60.6 m.矿体直接底板为含黏土质很高的 钙质黏土岩,厚 0.5 ～ 4.6 m;间接底板为薄层状泥质 白云岩,部分为角砾状泥质白云岩 ,矿层底板距奥陶 系下统 3.6 ～ 19.3 m.该区主要岩石力学参数见 表 1. 图 5 矿区地层柱状图 Fig.5 Geologicalprofileintheminearea 表 1 矿区主要矿岩力学参数 Table1 Mechanicalparametersofrocksinthemines 岩层 弹性模量 /MPa 泊松比 容重 /(kN·m-3) 内聚力 /MPa 内摩擦角度 /(°) 抗拉强度 /MPa 表土 50.0 0.30 18.0 0.08 20.0 0.05 第四纪沉积岩 800.0 0.23 21.0 0.32 21.5 0.26 灰岩 50 000.0 0.25 24.0 8.40 40.0 4.70 石膏 7 000.0 0.30 23.6 5.40 38.6 2.80 河北省邢台县尚旺庄石膏矿区的五个矿先后建 于 1984— 1998年期间,有第二石膏矿、太行石膏矿 、 林旺石膏矿 、邢燕石膏矿和康立石膏矿 (原恒昌石 膏矿)五座矿山 .塌陷区发生在康立矿和林旺矿之 间 ,塌陷区周边是各矿井的安全矿柱 (图 4).将塌 陷区所对应的采空区近似成矩形区域, 则坚硬顶板 长度为 2a=280 m, 宽度为 2b=180 m.矿柱群的总 横截面积占采空区面积的比率为 ζ=39.16%,矿柱 的高度 H=8 m.根据流变试验结果, 运用遗传算法 编程进行参数识别 [ 20--22] , 得到矿柱和顶板的流变模 型参数分别为 ηp1 =6.41 ×10 8 MPa·h, ηp2 =4.50 × 10 10 MPa· h, kp1 =27.7 ×10 3 MPa, kp2 =2.95 × 10 4 MPa;顶板流变参数为 ηr1 =1.03 ×10 16 MPa·h, ηr2 =4.68 ×10 12 MPa·h, kr1 =50 ×10 3 MPa, kr2 = · 146·