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·
第2期 王金安等：采空区坚硬顶板流变破断力学分析 ·147° 1.44X10MPa 采空区顶板黏弹性变形分析模型，通过弹性一黏弹 (1)作用在顶板上的总的均布载荷 性准则及LaPc数值逆变换，得到采空区顶板挠 49+Yh=21×145+24×38=3.957MP9 度随时间的变化关系. (2)顶板的抗弯刚度 (2)从考虑矿岩流变本构的矿柱顶板力学体 ER D 50×103×38 系分析可以看出，即使处于较为坚硬岩体环境的采 12(1-v)12X(1-0.25) 空区系统，在给定比较长的时间跨度，顶板位移仍会 228.63X10MPam: 达到相当大的量值.因此考虑岩体流变性质的稳定 (3)均布面力系数 性分析，不仅针对软岩采空区，在硬岩空区力学状态 4器39166 随时间变化的分析中也同样具有重要意义. (4)顶板破坏条件4，边界形成铰支时位移控 (3)由于影响矿柱顶板体系长期稳定性因素较 多，尤其是石膏矿，其中风化和水作用同样对采空区 制条件为 [oT]a 长期强度有着至关重要的作用，而本文忽略了这方 48D -=1.21 面的影响，因此应考虑风化，水作用下59对流变参 顶板中心破坏时Gm= 3x2D1 数进行折减，从而来反映这些因素的影响. %()= 参考文献 [o小，故有W=245四 【l】W ang JA Shang X C Ma H T nvesti讴ation of catastophic (5)将矿柱顶板黏弹性参数代入式(23，采用 giound collpse in X ingtai gypsum mines in China Int J Rock 数值Laplacej逆变换，得到固定，铰支下顶板中心点 Mech Mn Si200845(8):1480 位移图（图6，通过对图6分析，可以看出采空区 Shen NQ Yang JW ZhengX P Neural nework based collapse 顶板在大约6的时间短边破坏形成约11时顶 Predictin ofm ned out area CoalGeol Expkr 2001 20(3):42 板中心出现塑性铰. (慎乃齐杨建伟郑惜平基于神经网络的采空区塌陷预测. 煤田地质与勘探，200120(3为42) 45 4.0 [3 LaiX P Zhang L J CaiM F Application of neural nework o 3.5 铰支 the statistics and prediction of dynam ica l danage and evolution n 日3.0 the large sca le mine out area supported by ock based composite 25 顶板内部破断 materia ls JUniv SciTecnol Beijng 2003 25(4).301 20 固定 张1.5 (来兴平，张立杰蔡美峰.神经网络在大尺度采空区损伤演 1.0 边缘断裂形成塑性故 化统计与预测中应用.北京科技大学学报。200325(4 n.s 301) 1015 20 25 [4 W ang JA ShangXC LiuH etal Suudy on fmcuremechanim 时间.a and catastrophic oollapse of strong oof stmata above the m ned 图6固定、较支情况下顶板挠度曲线及破坏点 areg JChina Conl Sop 2008 33(8)850 F6 Deflec tin of the roofwith fxed and hinged boundares and (王金安，尚新春刘红，等.采空区坚硬顶板破断机理与灾 failure ponts 变塌陷研究.煤炭学报.200833(8,850) I5 BiswasK A unique appronch to de temine the time dependent n 应当指出，导致事故的直接原因是：尚汪庄矿区 siu stength ofcoalpillry/Proceed ngsof the Seond Intema tin 历经十多年开采，积累了大量未经处理的采空区，形 alWorkshop on CoalPillrMechanics and Desgn Vail 1999 5 成大面积顶板冒落的隐患：矿房超宽、超高开挖，导 AuvnyC Homand E HoxhaD The influence of relative humidit on the rate of convergence n an undeg ound gypsum mine ht J 致矿柱尺寸普遍偏小，稳定性较差；无序开采，在无 RockMech Min Sci 2008 45(8):1454 隔离矿柱的康立石膏矿和林旺石膏矿交界部位，形 [7 Caste llnza B GeropmatouE Nova R An atempt o predict the 成薄弱地带，受采动影响和蠕变作用遭到破坏，从而 failure tme of abandoned m ne Pillars Rock Mech Rock Eng 诱发了本次大面积采空区顶板冒落、地表塌陷事故. 200841(3,377 本文借鉴该案例，旨在为预测大范围采空区坚硬顶 I 8 StepiD GpdaG Viscoplastic behaviour around advanc ing tun 板变形和破坏提供新的分析思路和方法. nels n squeezing rocks RockMech Rock Eng 2009 42(2):319 I9 SivaniC De soyer T Bonelli Discre te malelling of tmede 5结论与讨论 pendent rockfill behaviour Int J Numer Anal Methals Geomech 200933(5665 (1)考虑岩体的流变性，建立了由矿柱支撑的 [10 Malan D F Smultng the tmedependent behavour of excava
第 2期 王金安等:采空区坚硬顶板流变破断力学分析 1.44 ×10 6 MPa. (1)作用在顶板上的总的均布载荷 q=q0 +γh=21 ×145 +24 ×38 =3.957 MPa; (2)顶板的抗弯刚度 D= Eh 3 12(1 -ν 2 ) = 50 ×10 3 ×38 3 12 ×(1 -0.25 2 ) = 228.63 ×10 6 MPa·m 3 ; (3)均布面力系数 = nA 4ab =39.16%; (4)顶板破坏条件 [ 4] , 边界形成铰支时位移控 制条件为 w0 ≥ [ σT] a 2 h 2 48D =1.21 cm, 顶板中心破坏时 σxmax = 3π 2D 2h 2 1 a 2 + ν b 2 w0 (t0 )= [ σT] ,故有 w0 =2.45 cm; (5)将矿柱顶板黏弹性参数代入式 (23), 采用 数值 Laplace逆变换,得到固定 、铰支下顶板中心点 位移图 (图 6).通过对图 6分析, 可以看出采空区 顶板在大约 6 a的时间短边破坏形成, 约 11 a时顶 板中心出现塑性铰. 图 6 固定、铰支情况下顶板挠度曲线及破坏点 Fig.6 Deflectionoftheroofwithfixedandhingedboundariesand failurepoints 应当指出,导致事故的直接原因是:尚汪庄矿区 历经十多年开采 ,积累了大量未经处理的采空区,形 成大面积顶板冒落的隐患;矿房超宽 、超高开挖 ,导 致矿柱尺寸普遍偏小, 稳定性较差 ;无序开采 , 在无 隔离矿柱的康立石膏矿和林旺石膏矿交界部位 ,形 成薄弱地带 ,受采动影响和蠕变作用遭到破坏 ,从而 诱发了本次大面积采空区顶板冒落、地表塌陷事故 . 本文借鉴该案例 ,旨在为预测大范围采空区坚硬顶 板变形和破坏提供新的分析思路和方法. 5 结论与讨论 (1)考虑岩体的流变性 ,建立了由矿柱支撑的 采空区顶板黏弹性变形分析模型, 通过弹性--黏弹 性准则及 Laplace数值逆变换 , 得到采空区顶板挠 度随时间的变化关系. (2)从考虑矿岩流变本构的矿柱--顶板力学体 系分析可以看出, 即使处于较为坚硬岩体环境的采 空区系统 ,在给定比较长的时间跨度,顶板位移仍会 达到相当大的量值 .因此考虑岩体流变性质的稳定 性分析,不仅针对软岩采空区 ,在硬岩空区力学状态 随时间变化的分析中也同样具有重要意义. (3)由于影响矿柱顶板体系长期稳定性因素较 多, 尤其是石膏矿, 其中风化和水作用同样对采空区 长期强度有着至关重要的作用, 而本文忽略了这方 面的影响 ,因此应考虑风化 ,水作用下 [ 5--6] 对流变参 数进行折减,从而来反映这些因素的影响. 参 考 文 献 [ 1] WangJA, ShangXC, MaH T.Investigationofcatastrophic groundcollapseinXingtaigypsumminesinChina.IntJRock MechMinSci, 2008, 45(8):1480 [ 2] ShenNQ, YangJW, ZhengXP.Neuralnetworkbasedcollapse predictionofminedoutarea.CoalGeolExplor, 2001, 20(3):42 (慎乃齐, 杨建伟, 郑惜平.基于神经网络的采空区塌陷预测 . 煤田地质与勘探, 2001, 20(3):42) [ 3] LaiXP, ZhangLJ, CaiM F.Applicationofneuralnetworkto thestatisticsandpredictionofdynamicaldamageandevolutionin thelargescalemine-outareasupportedbyrock-basedcomposite materials.JUnivSciTechnolBeijing, 2003, 25(4):301 (来兴平, 张立杰, 蔡美峰.神经网络在大尺度采空区损伤演 化统计与预测中应用.北京科技大学学报, 2003, 25(4): 301) [ 4] WangJA, ShangXC, LiuH, etal.Studyonfracturemechanism andcatastrophiccollapseofstrongroofstrataabovethemined area.JChinaCoalSoc, 2008, 33(8):850 (王金安, 尚新春, 刘红, 等.采空区坚硬顶板破断机理与灾 变塌陷研究.煤炭学报, 2008, 33(8):850) [ 5] BiswasK.Auniqueapproachtodeterminethetime-dependentinsitustrengthofcoalpillars//ProceedingsoftheSecondInternationalWorkshoponCoalPillarMechanicsandDesign.Vail, 1999:5 [ 6] AuvrayC, HomandF, HoxhaD.Theinfluenceofrelativehumidity ontherateofconvergenceinanundergroundgypsummine.IntJ RockMechMinSci, 2008, 45(8):1454 [ 7] CastellanzaR, GerolymatouE, NovaR.Anattempttopredictthe failuretimeofabandonedminepillars.RockMechRockEng, 2008, 41(3):377 [ 8] SterpiD, GiodaG.Visco-plasticbehaviouraroundadvancingtunnelsinsqueezingrock.RockMechRockEng, 2009, 42(2):319 [ 9] SilvaniC, DésoyerT, BonelliS.Discretemodellingoftime-dependentrockfillbehaviour.IntJNumerAnalMethodsGeomech, 2009, 33(5):665 [ 10] MalanDF.Simulatingthetime-dependentbehaviourofexcava- · 147·
。148 北京科技大学学报 第33卷 tins of hard ock RocMech Rock Eng 2002 34 (4):225 mericalmehods of inverse ngml transfomation in viscoelastic 【l川Napier JA↓MaknD E A visoplaste discontnuum malel of mechanics AdvMe吨199929(3为317 tie-dependent fmcuure and se ismic ity effect n brittle oock Int (魏培君，张双寅，吴永礼·黏弹性力学的对应原理及数值 RockMechMin Sci 1997 34 (7):1075 反演方法.力学进展199929(3片3引7) 12]FindkyW N LaiJS Ommn K Creep and Rekxa tion ofNan [19 PengF ChengY J Liu Y F et a]Numerical inversion of La lner Visooelastic Material with an Intioduction o Linear Vis Plce transfoms n viscoelastic prob kms by Fourier series expan coelasticit New York North Holknd Company 1976 s知Chin JTheor APPIM吨200840(2：215) 13]Cohen AM Num eric IMethods for Laphce Tmansfom Inversion (彭凡陈跃军，刘一凡等.基于Faurier吸数展开的L即ce New York Sprnger 2006 数值逆变换.力学学报.200840(2)：215 14]CaiME HeM C Lu D Y RockMechanics and Engineerng [20 WangX P Cao LM Theor of Genetic Algorith Aplication Be ijing S知ce Pres2002 and Sfware Realition X iap Xi'an JiaoongUniversity Press 蔡美峰，何满朝。刘东燕.岩石力学与工程.北京：科学出 2002 版社，2002) (王小平，曹立明.遗传算法理论，应用与软件实现.西安： 15]DulnerH Abate J Numerical nversion of Lapace transoms 西安交通大学出版社.2002) by re latng them o the finite Fourier cosine tansom J Assoc [21]Sandrne L G Expermenal Desgn Optmiation and Thema ComputMach 1968 15(1):115 Thysical Parameter Estma tion of Comn positeMa teria sUsing Genet [16]DuE Numerical nversicn oflapkce tmansoms an efficient c Agoritms[Dissermtion.Virgn Virgnia Poptechnic h mprovement oDubner&Abate smethal ComnputJ 1974 17 stime and S数e Universit议1999 (4):371 【2☑HeX Lis】LuY et a]Hentificatin of pemeabilit coeffi 17]Cnmp K Numerical nversion of Laplace transoms using a cientof rockmass n dam foundation based on genetic neutml net Fourier series appoxi ation J Asoc Comput Mach 1976 23 works Chin JRockMech Eng 2004 23(5)751 (1)片89 (何翔，李守巨刘迎曦等。基于遗传神经网络的坝基岩体渗 18]WeiP J ZhangSY WuYI,Corepondence prnc ple and nu 透系数识别.岩石与力学工程学报，200423(5,751)
北 京 科 技 大 学 学 报 第 33卷 tionsofhardrock.RockMechRockEng, 2002, 34 (4):225 [ 11] NapierJAL, MalanDF.Aviscoplasticdiscontinuummodelof time-dependentfractureandseismicityeffectinbrittlerock.IntJ RockMechMinSci, 1997, 34 (7):1075 [ 12] FindleyW N, LaiJS, OnaranK.CreepandRelaxationofNonlinearViscoelasticMaterials:withanIntroductiontoLinearViscoelasticity.NewYork:NorthHollandCompany, 1976 [ 13] CohenAM.NumericalMethodsforLaplaceTransformInversion. NewYork:Springer, 2006 [ 14] CaiMF, HeMC, LiuDY.RockMechanicsandEngineering. Beijing:SciencePress, 2002 (蔡美峰, 何满朝, 刘东燕.岩石力学与工程.北京:科学出 版社, 2002) [ 15] DubnerH, AbateJ.NumericalinversionofLaplacetransforms byrelatingthemtothefiniteFouriercosinetransform.JAssoc ComputMach, 1968, 15(1):115 [ 16] DurbinF.NumericalinversionofLaplacetransforms:anefficient improvementtoDubner＆Abate' smethod.ComputJ, 1974, 17 (4):371 [ 17] CrumpKS.NumericalinversionofLaplacetransformsusinga Fourierseriesapproximation.JAssocComputMach, 1976, 23 (1):89 [ 18] WeiPJ, ZhangSY, WuYL.Correspondenceprincipleandnumericalmethodsofinverseintegraltransformationinviscoelastic mechanics.AdvMech, 1999, 29(3):317 (魏培君, 张双寅, 吴永礼.黏弹性力学的对应原理及数值 反演方法 .力学进展, 1999, 29(3):317) [ 19] PengF, ChengYJ, LiuYF, etal.NumericalinversionofLaplacetransformsinviscoelasticproblemsbyFourierseriesexpansion.ChinJTheorApplMech, 2008, 40(2):215) (彭凡, 陈跃军, 刘一凡, 等.基于 Fourier级数展开的 Laplace 数值逆变换.力学学报, 2008, 40(2):215 [ 20] WangXP, CaoLM.TheoryofGeneticAlgorithm, Application andSoftwareRealization.Xian:XianJiaotongUniversityPress, 2002 (王小平, 曹立明.遗传算法理论, 应用与软件实现.西安: 西安交通大学出版社, 2002) [ 21] SandrineLG.ExperimentalDesignOptimizationandThermophysicalParameterEstimationofCompositeMaterialsUsingGeneticAlgorithms[ Dissertation] .Virginia:VirginiaPolytechnicInstituteandStateUniversity, 1999 [ 22] HeX, LiSJ, LiuYX, etal.Identificationofpermeabilitycoefficientofrockmassindamfoundationbasedongeneticneutralnetwork.ChinJRockMechEng, 2004, 23(5):751 (何翔, 李守巨 刘迎曦, 等.基于遗传神经网络的坝基岩体渗 透系数识别.岩石与力学工程学报, 2004, 23(5):751) · 148·