大型地下洞室地震灾变过程三维动力有限元模拟 (1)

第30卷第3期 岩石力学与工程学报 V01．30No．3 2011年3月 ChineseJournalofRockMechanic．YandEngineering March．2011 大型地下洞室地震灾变过程三维动力有限元模拟 张志国1～，肖 明1～，陈俊涛1t2 (1．武汉人学水资源与水电工程科学国家重点实验室，湖北武汉430072： 2．武汉大学水工岩石力学教育部重点实验室，湖北武汉430072) 摘要：地震作用是影响大型地下洞室群长期安全运行的r丰要外界因素。采用时程分析法模拟地震灾变中洞室围岩 地震响应过程，是研究地下洞室抗震的一种有效方法。为此，开发地下洞事地震灾变三维动力有限元数值模拟系 统。该系统采用显式中心差分法求解。为加快求解速度，对程序进行并行化处理并在积分方法上提出多高斯点分 区混合积分法；考虑高应变率下材料的强化特性和循环荷载作用下围岩损伤特性，提出适用于地下洞室抗震分析 的岩体动力本构模型；人工边界采用黏弹性边界；考虑地震波入射的方向性、多面性和非一致性，提出适用于地 下洞室抗震分析的地震波空间斜入射法。通过计算通用数值分析软件的验证算例，表明该系统程序求解的正确性。 同时，通过对鲁基场窑洞式地下厂房地震灾变时程模拟，表明系统程序在实际工程问题分析中的可靠性和实用性。 关键词t数值模拟；地下洞室；地震灾变；三维有限元法；动力本构模型；地震波输入 中圈分类号lO242 文献标识码lA 文章编号l 1000—6915(2011)03—0509—15 SIMULATloNOFEARTHQUAKEDISASTERPRoCESSoF LARGE-SCALEUNDERGRoUNDCAVERNSUSING THREE—DlMENSloNALDYNAMICFINITEELEMENTMETHoD ZHANGZhigu01～，XIAOMin91·2CHENJunta01·2 (1．StateKeyLaboratoryofWaterResourcesandHydropowerEngineeringScience，WuhanUniversity，Wuhan，Hubei430072， China：2．KeyLaboratoryofRockMechanicsinHydraulicStructuralEngineeringofMinistryofEducation，WuhanUniversity， Wuhan，Hubei430072，China) Abstract：Theearthquakeimpactisamajorexternalfactorinfluencingthelongtermstabilityofunderground caverns．Time-historyanalysismethodisabletosimulatetheearthquakeresponseofsurroundingrockin undergroundcavernsduringearthquakedisasterprocess．Ithasbeenprovedtobeaneffectivemethodtoanalyze theaseismicissueofundergroundcaverns．Thenumericalsimulationsystemusingthree—dimensionaldynamic finiteelementmethodisdevelopedtostudytheearthquakedisasterinundergroundcaverns．Centraldifference methodisemployedinthissystemtosolvetheproblem．Toenhancethesolvingspeed，parallelproceduresand hybridGausspointintegrationmethodareproposedinprogramming．Strengtheningfeaturesofmaterialsunder higllstrainrateanddamagefeaturesofsurroundingrockundercyclicloadingareconsidered．Thedynamic constitutivemodelforrockmassesispresented．Itissuitablefortheaseismicanalysisofundergroundcavems． Visco—elasticboundaryisadoptedinartificialboundary．Thespatialobliqueincidentmethodwhichissuitablefor theaseismicanalysisofundergroundcavernsisproposed．Itisabletoreflectthespecificincidentdirection，the multi—incidentsurfacesandtheinconsistencyofseismicwave．Anexampleisgiventoverifythecorrectnessofthe 收稿日期l 2010一II—15：修目日期l2010—12—26 基金项目t国家自然科学摹金莺人项H(90715042)：国家杰出青年科学基金项lq(50725931)：武汉大学博士研究生科研自事基金项目 作者誓介t张志国(1984一)，男。2007年毕业于武汉大学水利水电学院水工结构专业。现为博士研究生，主要从事地F洞窀围岩稳定数值计算方面的 研究上作。E-maih250566478@qq．colin 万方数据 ·510· 岩石力学与工程学报 2011年 developedsystembycomparingitscalculationresultswithcommonsoftware7scalculationresults．Acasestudyis thenconductedtosimulatetheearthquakedisasterprocessofthecave-typedundergroundcavernsatLujichang hydropowerplant．Theresultsderivedfromtime—historyanalysisindicatethereliabilityandpracticabilityofthe developedsystem． Keywords：numericalsimulation；undergroundcaverns：earthquakedisaster；three-dimensionalfiniteelement method；dynamicconstitutivemodel；seismicwaveinput ‘ 1引 言 随着国民经济建设的发展，我国西部水电开发 强度前所未有，已建、在建了一批大型、特大型水 电工程，均采用地下厂房的布置形式。我国西南部 是地震多发地，地震烈度大多在VII度以上。地震 作用是影响大型地下洞室群长期安全运行的主要外 界因素。研究地震荷载作用下地下洞室群围岩的损 伤、破坏、垮塌 机制 综治信访维稳工作机制反恐怖工作机制企业员工晋升机制公司员工晋升机制员工晋升机制图 ，建立大型地下洞室群地震灾 变过程数值模拟系统，对地下洞室的长期安全运行 及防止重大工程事故发生有藿大意义。 地F结构抗震研究涉及地震学、动力学、岩石 力学、喷锚支护理论等多项学科，一直是国内外的 研究热点。日本、美国等国外学者在地铁、城市地 下空间等浅埋土体工程抗震方面研究较多。我国基 金委将“重大工程动力灾变”作为重大研究计划后， 国内学者在已有研究基础上结合在建的大型水电工 程，在深埋岩体洞室抗震方面做了较多研究。如， J．B．Liu等【l’2j在地下结构动力时程分析的人工边界 及地震波荷载输入方法方面做了较多研究；李小军 等p“J在地下洞室群地震反应显式有限元求解方法 方面做了研究；钱七虎等15~71对岩石动力特性及本 构做了研究；李海波掣8】对地下洞室群地震荷载输 入方法做了研究；陈健云掣w采用相互作用法对溪 洛渡地下厂房震灾进行了分析；张志国等【IoJ对三维 弹塑性动力有限元隐式求解进行了研究。以往研究 多将重点放在人工边界及荷载输入、波动方程求解、 岩石动力本构等某学科领域内，分析平台多为商业 软件或在此基础上进行的二次开发，对地下洞室地震 灾变的系统过程模拟及自主程序开发研究较少。 本文着眼于水电站地下厂房洞室群地震灾变动 力过程模拟及系统程序开发，从工程应用的角度出 发，对地下洞室动态时程抗震分析方法进行了全面 研究，开发了地下洞室地震灾变过程三维动力有限 元数值模拟系统。本文程序采用显式中心差分法实 现对运动学基本微分方程的增量求解；采用OpenMP 并行模式对系统程序进行并行计算处理，并在积分 层面上提出多高斯点分区混合积分方法，在保障计 算精度的前提下，提升了计算速度；考虑高应变率 状态下围岩物理力学参数的提高和地震波动循环荷 载作用下围岩材料的疲劳损伤，提出了适用于地下 洞室抗震分析的岩体动力本构模型；波动方程求解 时，模犁边界采用人工黏弹性边界，并提出了适用 于地下洞室抗震分析的地震波空间斜入射法；通过 与FLAC3D数值分析软件的对比计算，表明本文程 序的正确性；针对不同实际地下工程抗震分析要求， 讨论了地震波的选取、滤波、人工合成及模型网格 划分原则；最后对鲁基场窑洞式地下厂房洞窜群进 行了地震灾变过程三维动态模拟。计算结果表明， 本文程序计算速度能满足工程分析要求，计算规律 正确，计算结果合理，可用于实际地下工程地震灾 变过程模拟。 2地震波动场的显式有限元求解 水电站地下厂房洞室深埋于山体之中，其地震 响应主要受地震波传播影响，属于近场波动问题。 工程区域内地震波动场的求解是地下洞室地震灾变 过程模拟的关键。对于复杂实际工程三维模拟，一 般采用有限元法实现对波动方程的积分求解。有限 元的求解方式，可分为隐式和显式2种。隐式方法 需组装分解总体刚度矩阵，受自由度限制对计算机 内存容量有较高要求。显式方法不需组装分解总体 刚度矩阵，对计算机内存要求较小。但其为条件稳 定算法，一般计算时步较小，计算耗时较长，计算 量大，对CPU有较高要求。从计算机硬件的发展看， 随着多核CPU并行计算的推广，计算速度成倍提 高，为大型地下洞室工程地震动力响应时程计算的 求解提供了坚实的硬件平台。从较成熟的岩土动力 分析软件看，如ABAQUS／EXPLICIT，FLAC3D， ANSYS／LSDYNA等，均采用显式解法。鉴于此， 本文程序采用显式方法完成工程区域波动场的求 解。 万方数据 第30卷第3期 张志国，等．大型地下洞室地震灾变过程三维动力彳r限元模拟 ·511· 2．1显式中心差分法 根据系统运动学基本定理，波动方程可表示为 增量拉格朗日格式的运动微分方程fll】： 删=／嘲一／“一，。娜 (1) 式中：石为网格结点加速度；fh，／嘲分别为网 格结点内力与外力；．厂蛳为阻尼力；M为质量矩 阵，在程序中采用集中质量矩阵形式，定义如下： M=Ip【妒]1【妒】d矿 (2) 式中：[纠为函数谚的矩阵，谚在分配给结点f的区 域内取1，在区域外取0。本文程序中，对于六面体 八结点 单元 初级会计实务单元训练题天津单元检测卷六年级下册数学单元教学设计框架单元教学设计的基本步骤主题单元教学设计 ，结点f的区域为其临近单元体积的1／8。 为便于矩阵求逆，本文程序应用陈育民和徐鼎平【l2J 所述的局部阻尼力的形式，将．厂4螂定义如下： ／‰p=口l厂矾一／hlsi嘶 式中：signa为网格结点速度方向符号；口为局部 阻尼系数，口=兀D，D为临界阻尼比，对于岩土材 料，D的范围一般为2％～5％。 对式(1)的显式求解方法有多种，本文程序采用 具有二阶精度的中心差分法(见图1)。 刚厂。，许～。 I．垒!：：!!：．1．垒!：：!!一l卜——————叫．_—=————．1 图1中心差分法时间轴 Fig．1Timelineofcentraldifferencemethod 假定时刻0，t1，f2，⋯，t”的位移、速度和 加速度均已知，其t”1时刻的位移、速度、加速度 可采用中心差分近似，即 a”+1=M。1(／嘲一厂“一／妇唧) (4) 西一牡=舀”一V2+Atii” (5) 口”+1=口一十△fn+l／2a”+1 72 (6) 其中， At=(At川72+Atn-112)／2，At”Ⅳ2=t“1一t“(7) 式中：口一1，a一分别为t”1和t“时刻的位移；矿， a”¨2，a一1分别为t一，t”1坨，t”1时刻的速度；西”， a”1分别为t一，tnt时刻的加速度。 显式中心差分法计算流程可描述为： (1)根据初始条件扩，扩计算：伊=M一(fo一 厂蛔)，a”¨2=西o。 (2)对每个时间步： ①由式(5)计算速度a”¨2； ②由式(6)计算t”1时刻的位移口”1； ③由t“1时刻的运动方程式(4)计算t”1时刻的 加速度西”1。 有限元显式计算是条件稳定算法，在大型地下 洞室地震响应模拟中，常常受到求解稳定和计算耗 时长的制约。这是因为要确保计算稳定，则一般需 要取较小的计算时间步长△f，从而增加系统循环步 数，延长计算时间，因此需要在确保计算稳定的前 提F，尽量选取较大的&。从波动有限元的物理意 义考虑，垃可理解为保障在一个计算时步内波的传 播距离不大于任何一个单元尺寸的最短时间。可采 用下式估算： 出≤口哑n苦(o．8。≤a一<0．98)(8) 式中：e为单元波速。对六面体八结点单元，，为单 元体积和单元6个表面的最大面积之比。 式(8)确定的系统时间步长，仅考虑地震波在计 算模型中的传播稳定性，并未考虑人工边界条件、 锚杆锚索等支护措施的计算稳定对系统时步的要 求。因此，在实际工程计算中，需针对不同情况对 系统时间步&进行调整。 为提高大型地下洞室三维模型的求解速度，缩 短计算耗时，提高程序的实用性，本文程序在求解 流程和积分方式上进行了优化。求解流程上，针对 显式有限元求解中各单元内力独立计算的特点，采 用可共享存储的、基于线程的OpenMP并行模式对 程序进行并行化处理，使本文程序在多核计算机和 大型工作站上的计算速度成倍提高。积分方式上， 针对显式计算中结点内力求解在每系统循环步中耗 时最多(约占80％)的特点，提出了单高斯点和八高 斯点分区域混合积分求解法。在显式有限元计算中， 单元结点力由各单元单独积分计算所得，单元之间 内力求解不相互干扰，这为单高斯点和八高斯点混 合积分提供了统一的数学格式。在实际工程有限元 模型求解时，根据计算任务对单元进行积分区域划 分，对重点分析的计算区域、边界单元等部位采用 八高斯点积分，确保计算精度及人工边界的数值稳 定。对其余模型区域采用单高斯点积分，并进行沙 漏修正。该处理方法在保障计算精度的前提下，简 化了计算流程，提高了计算效率。2种高斯点分区 积分会对结点内力计算产生一定误差，从而影响整 万方数据 ·512· 岩石力学与工程学报 2011芷 个系统的计算稳定，在理论上看需适当减小系统时 间步长。从本文程序对溪洛渡、瀑布沟、渔子溪、 鲁基场等众多实际大型地下洞室群的地震响应模拟 分析情况看，尚未出现因此而需减小系统时步的情 况。 2．2结点内力求解 在式(1)的显式求解中，最关键的是结点内力 ．厂觚的求解。在增量拉格朗日显式计算中，首先采 用率本构方程由系统结点速度求解计算时步内单元 应力增量，并对单元进行应力更新。然后将更新后 的单元应力积分到单元结点上，从而求得结点内力。 (1)变形描述 在增量拉格朗日显式计算中，变形采用应变率 ％描述： 白=毒=三(考一等]+三(誓+鲁]_ao．+‰ (9) 式中：纯为旋率张量；见为变形率张量；M为速 度分量；z，为坐标分量；f；，=l，2，3。 (2)应力更新 积分率本构方程的数值算法称为应力更新算 法。时刻t+dt的应力％O+dt)可以通过对应力率 线，积分得到： ％(f+dr)=％(f)+吃dt(10a) 柯西应力率可表示为 岛=醇+％铱+％q (10b) 式中：茚为焦曼应力率，醪=％巩，％为弹 性本构张量。 将式(tOa)求得的应力％(f+dt)积分至单元结 点，得到结点内力矢量厂h。 3 围岩动力弹塑性损伤本构 3．1洞室围岩动力弹塑性损伤本构模型 洞室围岩的动力损伤本构模型是地下洞室地震 灾变过程模拟的基础。这需要在理论研究的基础上， 结合大量的室内、室外试验，建立科学实用的岩石 动力损伤本构模型。由于这项工作的复杂性，目前 国内外尚未有普遍认可的适用于地下洞室围岩地震 分析的本构模型。鉴于此，笔者基于对地震灾变中 地下洞室围岩特性认识，在Mohr-Coulomb屈服准 则基础上，提出一种适用于地下工程围岩地震灾变 分析的本构模型。 与准静态洞室开挖过程相比，地震灾变中围岩 承受地震波循环荷载作用，表现出围岩材料的动力 强化特性和疲劳损伤特性。研究结果【5H副表明：动 荷载作用下，围岩应变率提高，对围岩材料的弹性 模量和抗破坏性能均有所提高；循环加载情况下， 围岩疲劳损伤，对围岩材料的弹性模量和抗破坏特 性均有所降低。两因素相互对立，又同时存在，在 本构模型中均应考虑。 (1)对于岩石动力强化特性的研究主要分为高应 变率(毒>l02s-1)和中应变率(10_4s-1<童<10-1$-1) 研究。其中，高应变率主要针对冲击爆破情况，由 于加载方式简单试验相对较多。而地下洞室地震作 用应变率一般在10一s-1<营<10—1s_。，属中等应 变速率，试验较少。从已有的研究成果【5·14～17】可以 得到以下结论：①岩石动力强度及弹性模量随应变 率的增加而增加；②岩石动力强度随应变率增加， 主要是因为黏聚力随应变率增加，而内摩擦角不受 应变率影响；③中等应变率的岩石动力特性受加 载方式、加载速率、岩石类型及试验设备影响，各 种试验结果存在较大的离散性，尚难在试验基础上 统计出一致规律。 (2)岩石损伤本构研究相对较多，其中损伤描 述主要有一、二、四阶等张量描述形式。针对地震 波传播的有向性，本文程序采用二阶损伤张量描述 地震中围岩受损情况。 考虑地震中围岩上述2种特性，本文将岩石动 力本构模型围岩动弹性模量表示为 最=P(g)E (11) 式中：E为围岩静弹性模量；P(舌)为大于1的应变 率函数，根据试验结果确定。 屈服准则采用带拉伸截止限的Mohr-Coulomb 屈服准则(见图2)，且假定压应力为负，拉应力为正。 线4B为剪切屈服准则f。=0，线BC为拉伸屈服准 则厂‘=0，即 厂s2I一鼍％+2旭p弦√％}(12a) 厂t=吠一∥ j 其中， 以=等 (12b) 万方数据 第30卷第3期 张志国，等．大型地下洞室地震灾变过程三维动力有限元模拟 ·513· 图2 Mohr-Coulomb屈服准则 Fig．2Mohr-ColIIombyieldcriterion 式中：C，缈分别为黏聚力和内摩擦角；西，《分 别为损伤后的有效第一、第三主应力；仃t为极限抗 拉强度，其最大值叱=cltan口o；T为损伤影响系 数，二阶张量形式下，T=1-40：+厦+研，D。， D2，D，为各主应力方向损伤系数，由损伤演化方程 求解。 每计算时步内，L Q(g)为常数，按照围岩非 关联性，屈服准则相应的势函数为 ，j一《以} (13a) gt：《 I 、 ’ 其中， 以=面l+sinv (13b) 式中：沙为围岩的剪胀角。 计算中单元应力采用考虑损伤后的有效应力形 0 0 1 I 吒(1一砬向)01 0 o-3(1一D3h)j (14) 受拉情况下取h=1．0，受压情况下取h=0．2t¨J。 围岩的损伤破坏过程是伴随其微裂隙发展和微观空 隙率增加，材料强度降低的过程。从宏观上看，围 岩损伤主要是由于受张拉而出现开裂破坏造成的。 根据一些实际工程的观察结果表明，大部分围岩的张 性破坏是围岩实际应变超过极限拉应变所致ll州。针对 三维情况，可以认为在单元应力进入塑性，且第一 主拉应变超过极限拉应变，则发生损伤破坏，即 f 0 D={． 【办(勺) 篙薹罱凳瓮 ∞，(E≥【￡】和／5>0) 、7 式中：矗(岛)为三维损伤演化方程，毛为围岩第一 主拉应变，【￡】为围岩极限拉应变。 对于在实际工程中未测试围岩的极限拉应变的 情况，在本文程序中用岩石抗拉强度R和弹性模量 E的比值代替，即 M=足／(kE) (16) 式中：k为安全系数。 在损伤区，围岩应力和物理参数都随累积的塑 性偏应变的增大而减小，累积塑性变形越大，围岩 损伤程度越高。在围岩材料强度极限前的相邻区域 内微裂纹损伤增加是非常迅速的，对于这种损伤破 坏的变化速率通常可以用指数函数来描述。因此， 三维损伤演化方程可以用下式表示： ，r———一、 B=l～expl—R√掣掣J(f=1，2，3)(17) 式中：卵为第f个主应变方向上的塑性偏应变，R 为材料的损伤常数。 为便于在后处理中显示围岩损伤状态，在结果 文件输出时，将二阶损伤张量转换为一阶损伤标量， 表示为：D=、／胼+历+碍。 3．2应力修正及损伤计算 单元应力计算是显式有限元计算的核心。在 m+l时刻弹颦性损伤计算中，首先根据节2．2所述 采用单元材料初始参数计算各高斯点弹性应力，并 判断单元是否处于损伤状态∞泐)≠0，此时单元损 伤系数为tn时刻计算所得)。若单元处于损伤状态， 则修正各高斯点应力为有效应力。在此基础上，考 虑损伤对屈服面的影响，进行屈服判断。若发生屈 服，则进行塑性应力修正，并且判断第一主拉应变 是否超过极限拉应变。若超过，则认为单元发生损 伤，并计算损伤系数。 单元塑性应力修正有多种方法，本文程序采用 使用手册120J中所述的Mohr-Coulomb屈服准则塑性 应力修正方法。该应力修正方法，适用于屈服函数 八吒)是线性函数的情况。根据塑性应力需在屈服 面上流动的原则，应保障tn+l时刻的塑性应力有 f(o"n+△吒)=0 (18) 对于线性函数，式(18)可转换为 厂(吒)+厂(△吒)=0 (19a) 厂(‘)=厂(·)一f(0。) (19b) 式中：厂(·)为屈服函数减去其常数项／(0。)。 q^— O O 0嘎 P．．．．．．．。．．．．．。。．．．。L =1J盯 式 万方数据 ·514· 岩石力学与工程学报 基于该变换，使用手册【20】推导了线性屈服函数 塑性应力修正公式： ∥井鸽(鲁](20， 式中：名为塑性流动因子；S为弹性本构：叫，掣 分别为弹性假定应力和塑性修正后应力。 本文程序本构模型在每计算时步f内，Q(g)为 常数，屈服函数为线性函数，可采用式(20)对考虑 损伤后的单元有效应力进行塑性修正。对于塑性剪 切破坏，有 叫=酬’一力5(％一钙Ⅳ0)l 叫=一4一旯8％(1一Ⅳy)} (21a) d=叫‘-)t8(一％～+呸)J 名s： ￡熊：：垡：! (2lb) (％一a=N矿)一(一％～+az)N矿 、 ’ 其中， q：K+芸G，a2：K一罢G(21c) 式中：K，G分别为体积模量和剪切模量；《+为损 伤修正后的等效弹性假定应力。 对于极限拉破坏，有 带叫刊+可’)鲁 秽科刊‘可‘)等 秽=∥‘ (22) =P--P1-矿exp五(笔)卜B。州，=一R厄丽)l 岛=(础州)+咣。+1)+钱州))／3(24) 3．3数值计算流程 已知所时刻模型结点速度矢量为v(tn)，系统 节点更新后的坐标为x(tn)，单元应力为cr(tn)，则 tn+l时刻系统计算流程如图3所示。 多 核 并 行 计 算 图3船H时刻计算流程图 Fig．3Flowchartattn+ltime 4人工边界 8 高 斯 点 循 环 计 算 人工边界实质上是工程区域内波动场有限元求 解时，为实现无限介质的有限化处理而人为虚设的 计算边界。人工边界是近场波动数值模拟的关键， 20世纪60年代以来，国内外学者已经进行了广泛 而深入的研究。国内的廖振鹏等【2卜231对局部人工边 界，尤其是黏弹性人工边界的发展做出了较大贡献。 刘晶波等弘2J已对黏弹性人工边界理论及物理模型 做了详细论述。以下主要讨论显式有限元计算中黏 弹性人工边界的计算格式及相应地下洞室地震波荷 载的输入方法。 4．1黏弹性人工边界显式格式 黏弹性人工边界为应力型人工边界，可以采用 在边界结点3个自由度上各施加一组弹簧+阻尼器 系统模拟。模型边界结点受力根据结点运动状态， 按下式求解： (1)法向方向，有 万方数据 第30卷第3期 张志国，等．大型地下洞室地震灾变过程三维动力有限元模拟 ·515· Pn=-'fan=彳(KN“+CN西) (2)切向方向，有 辟=彳听=彳(KTu+CT西) 其中， ，，，n 室群一般位于深埋岩体中，在地震灾变模拟中，不 ‘。 能像分析地表建筑物(桥梁、大坝等)一样，只考虑 地震波从底部入射，而应考虑从模型多个迎波边界 (26) 入射(见图4)。 疋，：比．旦．k：缘旦1 “ “R’ 1 1 R} (27) 、7 CN=pCp，cT=pcsJ 式中：KN，矸分别为法向和切向弹簧刚度；aN，aT 分别为法向和切向修正系数；G，cT分别为法向和 切向阻尼系数：U，zi分别为边界结点位移、速度； cP，C分别为介质压缩波与剪切波波速；P为介质 密度：冗为波源至人工边界点的距离；彳为边界点 的代表面积。 本文程序质量矩阵为集中质量矩阵，各网格结 点代表其周围真实质量的分布。为确保人工边界计 算中采用的结点代表面积与结点代表质量的一致 性，对于八结点六面体单元，边界点代表面积彳为 与结点相关的边界网格面面积的1／4。 对于入射边界，根据波场分解原理加时刻边界 点的运动状态是入射波场与散射波场的叠加。其中， 入射波场可根据入射地震波时程曲线直接计算求 得。散射波场为自模型内部的外行波，由模型内有 限元时步计算求得，在中心差分法中可认为是上一 时步入射边界点的运动状态。对于透射边界，根据 波场分解原理仂时刻边界点的运动状态是自由场与 散射波场的叠加。其中，散射波场与入射边界的散射 波场求解相同，自由场可采用傅氏变换计算求得。 从式(19)可看出，黏弹性人工边界通过在边界 点上施加集中力，引入到式(1)的求解中。因此，可 将人工边界作用力看作是施加在模型边界的外力， 从而简化计算，增强程序计算模块间的独立性。但 是，黏弹性人工边界计算在时域上的离散与中心差 分法在时域上的离散相差半个时步，从而引入了潜 在的数值稳定问题。这种时间离散上的不一致丰要 体现在，采用仂时刻边界点上的外行波场计算出的 边界力，作为边界荷载参与m+l时刻的计算。采用 较小的计算时步，可以避免边界计算数值不稳定情 况的发生。另外，通过对众多实际工程的计算看， 对边界单元进行系统时步的估算时，式(8)中口值应 小于0．5。 4．2地震波荷载输入 地震波的荷载输入方法主要是研究确定节4．1 中所述入射边界上的波形函数。水电站厂房地下洞 图4地下洞室地震波入射示意图 Fig．4Sketchofincidentseismicwavesinunderground CaVern 在图4中，由于地震波到达OA，OB两入射边 界上各模型网格结点的时间不同，因此在同一时刻 入射边界上各结点波动状态不同，在地震波输入时 应考虑这种不～致性。 在地下洞室抗震分析中，为考虑地震波的多边 界斜入射和边界波动场的非一致性，本文提出了地 震波空间斜入射法。该方法认为，地下洞室模型分 析区域相对整个实际地震区域是微小的，可以看作 质点，工程区域的地震波可以看作平面波。进而， 可以认为与地震波传播方向垂直的平面(图4中PQ 平面1内各质点波动状态是一致的。 定义：地震平面波最先到达的模型边界点为输 入基准点(图4中点D)；与地震波传播方向垂直， 且通过输入基准点的空间平面为输入基准面(图4 中PQ面)。 假定：根据胡进军和谢礼立【24J的研究成果，可 假设地震波在较小的工程区域内传播时，仅地震波 幅值受岩体阻尼影响，呈线性衰减。 设基准面地震波时程为∽=彬(f)(净1～3)，则 如图4中0．4，OB入射边界上各结点的波动时程为 u=彬(f一垃)(1一肛’} (28) At=L／C J 一 式中：&为地震波从基准面到边界入射结点的传播 时间，￡为入射结点到基准面距离，C为地震波波 速，∥为地震波幅值沿传播距离衰减系数。 在地震时程计算中，首先需要根据地表监测地 震波，反演输入基准面地震波及地震波衰减系数 ∥；在此基础上，采用式(20)计算模型输入边界结 万方数据 ·516· 岩石力学与工程学报 2011年 点地震波形。结合图4，具体过程如下： (1)设地表监测点为点D。建立整个工程区域 包括地表监测点、断层、不同岩层等在内的大模型 OEFB。该大模型网格主要用于对区域波动场的反 演，网格尺寸可设置较大，但需满足地震波频谱分 析的计算频域要求。 (2)将工程区域选取或合成的地震波 Ui=彬(f)作为监测点D的地震波。则大模型输入 基准面尸Q上的地震波为Ui=形(f)(卜鹏)，进一 步采用式(20)计算模型输入边界OB，OE上各结点 地震波。 (3)以点D已知地震波形为收敛控制 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 ，对 大模型进行弹性有限元波动求解，反演地震波衰减 系数p。 (4)建立用于洞室弹塑性损伤动力时程分析的 小模型OACB。根据反演所得地震波衰减系数∥及 基准面波形，计算小模型入射边界结点波形，进行 三维弹塑性损伤动力时程计算。 上述地震波荷载输入法考虑了地下洞室结构地 震波的入射特性，反映了地震波的多面入射和边界 波形的非一致性特点，使地下洞室抗震模拟更为真 实。笔者对空间斜入射与底部入射做了多项数值对 比试验，由于篇幅限制不再详述。试验结果表明， 地下洞室采用空间斜入射法较底部入射法，模型破 坏区、损伤程度会有明显增加。为保障工程安全， 避免重大事故，从偏安全的角度考虑，笔者建议在 地下洞室地震灾变仿真中，地震波输入应考虑入射 多面和非一致性。 5地震波荷载 5．1地震波的选取 科学的地震波选取是地下洞室群地震灾变过程 模拟的前提，直接影响计算结果的可靠性。目前， 对于水电站地下洞室群地震波的选取或人工合成并 无 规范 编程规范下载gsp规范下载钢格栅规范下载警徽规范下载建设厅规范下载 说明，学术界亦未有普遍的统一认识。在实 际工程分析中，常按照规范125】中方法选择或人工合 成地震波。虽然该规范主要针对地面建筑抗震 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 ， 但其在地震波选取中的基本原则与思路，值得地下 洞室在抗震分析时借鉴。 (1)由于未来地震的不可预测性和实测地震波 的多样性，可根据有限样本容量下的计算结果来估 计地震效应。水电站地下洞室一般位于崇山峻岭之 中，场地历史实测地震波一般较少。在具体工程分 析时，应全面分析工程区域内的震源构造，搜集相 关震源历史地震记录资料及实测地震波。同时，根 据地形地貌、岩层岩性、断层分布等，类比选取具 有相似场地的国内外实测地震波。这些地震波组成 初始大样本。 (2)根据地震动三要素：频谱特性、有效峰值、 持续时间，选取对地下洞室影响较大的几组地震波 组成小样本。其频谱特性应满足场地类别与设计地 震动分组要求，同时尽量选择对地下洞室影响较大 的低频波。一般来说，有效峰值和持续时间越大， 地下洞室的地震效应越大。小样本容量一般不小于 3组。 在实际水电站地下洞室群地震灾变模拟中，由 于其三维空间地质构造、洞室结构等分布的非对称 性，一般应选择输入x，Y，z三向地震波。同时， 需要说明的是，地震波的选取与人工合成对分析人 员的工程经验有较高要求，实际工程分析时在条件 允许情况下，建议采用专业地震工作者或地震局提 供的用于特定工程的实测或人工合成地震波。 5．2波形处理与网格尺寸 计算输入地震波的合理性直接影响计算结果的 精确性。因此，在计算前应对地震波进行基线校正 和滤波处理。基线校正是为了消除实测地震波计算 中出现的模型飘逸现象；滤波处理是为了消除数值 计算中地震波的无效频带，保留有效频带，确保计 算精度。地下洞室地震波的有效频带确定原则： (1)若显式计算中每时间步距为&，计算持续 时间为瓦=NAt的输入数字荷载的理论有效频带是 1／Tg--"1／2At【4】o (2)为保障离散时间步垃能反映地震波形，应 满足at／K≤o．1(其中，瓦为输入地震波中的最短周 期t41)。 (3)对输入地震波形进行功率谱分析，保留包 含90％能量的频段。 (4)输入地震波频段应满足模型网格的精确分 析要求。有限元网格模型建成后，其能分析的有效 频段即已确定。一般情况下可认为，其模型有效频 r r 段为而每～而每(矿，驴分别为模型最大、 最小单元特征尺寸)。对于本文程序采用的黏弹性 人工边界，若输入频带超出模型有效频段，则可能 出现“低频漂移”和“高频振荡”现象，影响计算 精度。 万方数据 镕30卷第3期 张志田．等大型地下祸室地震灾变过程=雏动力有限i模拟 袖A地震波频率与模型网格尺寸有密切关系。 两者不但存在上述(4)中的频段关系．为确保计算稳 定，网格尺寸太小从物理意义上讲，应大于波在一 个时间步距的传播距离，这也是式(5)的确立原则。 为保障计算精度．规范1261规定，模型网格在波传播 方向的尺寸应在所考虑最短波长的1／12～l腮范围 内取值。笔者曾就该尺寸范围进行了大量数值模型 试验，表明碍格尺寸上限难以超越1月。 6程序验证 为说明本文程序计算的正确性，本文分剐针对 程序的本构模型和人工边界进行验证。 6．1弹皇性掘仿奉黼证 奉文程序所用届服准则是在Mohr-Coulornb屈 服准则基础上，考虑了围岩动力强化和疲劳损伤软 化特性，本质是对Mohr-Coulomb屈服准则的一种 修正。为验证计算程序的正确性，采用本立程序对 使用手册脚1中FLAC3D所采用的薄板圆孔验证算例 进行计算，井将计算结果与FLAP”计算结粜进行 比较。需要说明的是：本文程序虽主要用于求解动 力时程问题，但通过增大局部阻尼，采用动态松弛 算法也可求解静力问题；为保障与FLAC“中Mohr- Coulomb本构模型计算条件的一致性，在该算例中 本构模型分别采用考虑疲劳损伤特性和不考虑疲劳 损伤特性2种。 算例l：腔设无限域薄板中赋存三向韧始应力 均为30Mp“压为负)，内有一圆孔开挖，计算孔 周应力分布。计算模型见图s(柚。材料参数见表1。 根据计算采用的奉构模型，可分为：工况l，本文 程序不考虑掘伤的本构模型(P=l，Q=l-R=0)： 鲨 二8 『*。 目 ‘d)m*§&*≈ 目5算例l计算结果 Fig5 Cal州lating惜u】Bof～ple 表1薄板圊孔算倒材料参数 Table1 P一㈣ofmartialin deofthinplatewith⋯k”hole 工况2，本文程序考虑损伤的奉构模型(尸=l，口=1- R=1k上况3，FLAC”的Mdar-Coulorab本构模型。 “)从工况】的计算结果看，塑性区分布范目 及应力分布规律均较合理(见圈5嘞，0))。第一主 应力照人值位于塑性区与弹性区的交界处。 №l!¨-¨_ 万方数据 ·518· 岩石力学与工程学报 (2)3种工况计算后的单元相对应力a．／Po，叫 R(R=一30MPa)沿距孔壁距离分布情况如图6所 示，工况1与3基本相同，差值在1．5％之内。工况 2第一主应力在峰值点前较工况l小，主要是因为 模型考虑损伤后，塑性破坏区应力部分释放。 芒 b R 搓 霞 垩 距孔壁距离，m 图6相对应力与距孔壁距离关系曲线 Fig．6Relationcurvesbetweenrelativestressanddistance fromholewall (3)从图5(d)可以看出，由内向外损伤系数逐渐 减小，规律合理。 综上，通过对3种工况的对比计算分析，说明 本文程序本构模翠是合理的，计算结果是正确的。 6．2人工边界验证 为验证本文程序黏弹性人工边界的正确性，同 样采用使用手册【201中FLAC3D的验证算例进行计算。 算例2：一根水平的弹性杆，长50m，宽lm， 弹性模量为25．7MPa，泊松比为0．286，密度l000 kg／m3。杆的左端施加一竖直向冲击荷载，右端为自 由面，左端施加人工边界。计算不考虑初始应力、 自重等。具体计算分3种工况：工况l，黏性边界 不考虑阻尼懈=0)；工况2，黏性边界考虑阻尼懈= O)；工况3，黏弹性边界不考虑阻尼似≠0)。算例2 模型见图7。 l m／s●．人I：边界 自In而I／＼r[1lllllIlllllllllllllllllllllllllIIIlItIIIII[1111I、，1 DL’、底部监测点 巾氤监测点 琐部监测点7 图7算例2模型 Fig．7Modelofexample2 记录分析杆底部、中心、顶部3个典型结点的 速度响应(见图8)。工况1，不考虑阻尼时，底部点 的速度响应与输入荷载完全一致。各监测点波形相 ， }∞ ● 吕 ￥ 2．0 1．5 ’∞1．0 喜0。5 O．0 一O．5 2．5 2．O ，I．5 · 1．0 暑 ≤O．5 O．0 —0．5 (a)工况1 (b)工况2 ，／s (c)工况3 图8算例2监测点速度时程图 Fig．8Velocitytime-historiesofmonitoringpointsofexample2 差时步与波的传播时间相同。顶部点自由面速度幅 值是输入荷载的2倍。最后2个脉冲是从模趔自由 面反射回来的波，反射波传播至底部后，模型基本 结束振动。计算结果的规律、量值与FLAC3D基本 完全相同。工况2，考虑阻尼时，波的传播规律在 工况1的基础上，表现出明显的衰减性。工况3， 采用黏弹性边界后，波的传播规律与工况l基本相 同，但在反射波传播至底部后，模型呈现整体向下 运动的趋势。这主要是因为，输入脉冲荷载向上， 模型整体向上移动，脉冲作用完毕后，模型在弹性 边界的作用下，恢复初始位置，向下运动。 综上可见，本文程序人工边界设置是合理有效 的。 7工程实例 7．1工程概况 为说明本文程序的实用性，选取位于云南普渡 河干流上的鲁基场窑洞式地下厂房进行实际工程分 析。该地下厂房洞室群由主厂房、尾闸室、引水洞、 母线洞、尾水洞等洞室组成。共3台机，从窑洞口 k 点，Zm如 部 、、一I●， 顶“， 一^挎=====一涨＼===j。遮 5 O 5 O 5 O 5 2 2 l l O O O 万方数据 蒉；30卷第3舰 张忠目．萼人型地T雌地《女皇塾里三丝生生塑里墨堡丝 『；u内依次分靠。受J区地彤地质条件嫩制．主厂房 _f『i置于大坝左岸下游【【J体内．为窑洞式。工程区地 震烈度为VIII度．存在地卜J房和高边坡抗震问题。 土厂房最大高度3725in．最大跨度24／n。洞 窀群处j。寒武系r统(∈1)弱风化白云质获岩中．岩 质。旺硬，裂隙中等发育。匿岩主要为¨l娄，局部 lv娄固岩。物理力学参数见表2。本算例．暂不考 虑镐杆、锚索、衬砌等支护措旖。由于本工程末进 {T现场地应力测试，本文计算初始地应力场采用考 虑河粹及开挖边坡影响的自晕应力场。围岩阻尼采 用局部阻尼的形式，取0157I。本上程未进行岩体 变形模量和黏聚力随应变率变化的试验．从工科安 全角度考虑．本算例尸’0函数均取I。 &2冉摹场地TJ房IH*物Ⅻ力学参数 Table2 PhysI⋯ch¨1c“MmcⅫ，。frockmsⅫof undomrO帅d⋯1 考虑到地下』房距地表较近，工程区域地震波 动场反演模掣和地下厂房洞室分析模型可采用同～ 模型。建立宵限元计算模型包括土厂房、尾闸窜、 尾水洞、母线桐及引水洞．共272304个六面体八 节点单兀，包古301l】2个结点，主厂房轴线A坩 为P冉向．』轴‘，p轴睡直指向F游，z轴与人地坐 标煎合。沿X．Y，：轴3个方向的计算范崮分别为 1760，1090．1410m。鲁摹场地F厂脐删窜群见 留9。为反映三维输入特性．奉文地震波采用汶川 地震中卧龙台记录的三向地震加速度时稚．井截取 毖中20s进行计算。按照节5中所述地震波波形处 理及网格尺寸确定原则，本工程滤波频率范围0l～ 100Hz．相应最大网格特征尺寸不大于7．0m。对 II|_繁 —删姗¨ k17自№aR日a 月10波⋯地震目龙台宴测加速度时程曲线 Fi910Measuredtime-hIm。wcⅢⅧ“一I∞【i舯 inWoI硼Bstati∞oh舯酣“ngWenchu“ cmqLlake 7．2工曩匠嚏谴动墙反演 根据节4,2所述波动场反漓疗浊．认为上述处 理后的卧龙波为丹摹场地Fr晴上程区地表监测点 实测波彤．作为反演试算收敛拄制标准，采用本文 程序弹性车构．计算反演地震波幅值随距离的衰减 万方数据 §6力々bI程学报 20ll《 系数口腱输入丛准面地震波形，需要说明的是．卧 儿强震台附近不存在大型地下洞室，实测地震渡来 受#t体tp涧室空洞影响．因此笔者认为在反{寅地震 波幅值鼬距离的衰减系数口及基准面输入地震被 形时．所用模型鹿为洞宣开挖前模型．地表监测点 反演及实测地震波位移时程曲线见圈11．反演壬丑关 系散092。波形吻台较好。地表监测点距输入基准 面1700m，考虑地袭自由而的放大作用．褒减系 数占=00279， 驴朋八|；||L．型．；y￥这 蹦II№＆＆■^&渣&女■地l波俺g时捏自线 F慵II■n1}hi鲫cl州8吖n口酬^T-dm劬b删 d1∞hⅫm临m辨Lmd川帅ngpoint "地下一童地l虫壹过曼■挂 (1)嗣宣丌挖模拟 在对地下洞率进行弹颦性损伤地震宄变模拟 前．首先需螫模拟整个地FJ房州室群开挖后围岩 应力扰动、塑性开裂艘坏、围岩拭伤婶Ilj岩状态的 扰动。本文程序通i：上调整局部#H尼，采用动卷松弛 法IU以求解静力问堰。奉，妾例中．请艘挫f}局部阻 尼为n8．对幽9所示洞审扦挖模型进行次性无支 护开挖。开挖后．洞周围岩塑性擞伤I曩圩嘶见圈12。 ：。。。．； ：i目 =!!I；：品I 守：卜 目12洞％”挖3’机组霞阿#避性捌伤E分布 Fi912DiStributionofpl*ticdamagezonc“unit啦6咖柠 an“exca％1lionslmulaIIon 涸审静山丹挖完毕后，土J脐两侧高边墙围岩 部分进入塑性损协区，深度住5m2内．按伤系数 城大达569x10～。母线洞、Jt闸宣等嗣事较小． 围岩破坏区较小。总体打泉．厂房开挖后洞周围岩 相对稳定．需璺说明的址．本文湖性损伤区蔗指根 据弹塑性授伤本构求解计算所得的，麻力状态进入 塑性且发生损伤破坏的区域。 (2)洞章地震灾变模拟 在洞审开挖围岩损伤、戍力扰动的摹础上，采 用节7．2反演的基准面地震波及动力衰减系数卢， 计算模拟地下洞室群地震灾变址程。术模型按约27 万个A呵体八结点单元．动力计算持时25玎前20s 有地震波输八1．采用配簧Intel(R)Core(TM)i5280 GHz4拨CPU的营通办公电脑计算耗时约108h． 计算速度能满足工程分析要求． 计算结果分析如下： ①从圈13可以肴出．各监测点位移时程规律 与反演输入地震渡基本丰fI同。由下临空面效应．各 临测点位移峰值均较输入地震渡肯大幅增加．其中 主J‘虏项拱监测点A增大80％。上游边墙监测点e 增大55％．下游边墙监测点D增人45％，曲坡监测 点F增大58％。临空面的放大效应．对表层田岩稳 定不利。在25s持时计算完毕时．随着地震波的外 行和阻佗耗敞，模掣系统恪体动能蕈木为0，选到 稳迎状态。 瑚13动力计算过程监铡点位移时艘曲线 Fig13nm}h㈣∞ofdisptar．cm∞tsofmonitoring pointsduringdyrmmieealeuhtionInoeess 在波形时间差、空间轮廓、塑性损伤及阻尼耗 散的影响下．各监测点位移时程曲线在鬻值上冉差 别．出现相对运动。这也是洞室在地震作儿IF失稳 的一个重要原因．由豳13可见。丰厂房}游边墙与 下辩边墙测点最大相对位移达6gm，约为输入渡肜 荷载幅值50％。考虑厂房帆组安装回堪泥稚上，能 在一定程度上减缓这种情况。 @从图14可以看出，在地震波0i我的作用下· 圉岩第三主应力处于剧烈的椭环震筠中．这种循环 壤-5苦-叶衅_毫JJ。 万方数据 第30卷第3M 张志田．等大攫地F相窀地震灾变过程三维动力冉限元模拟 主0．1谚譬i ¨● 棚f*∞* 目14目周围岩典型∞位第!±成力时程∞拽 F嘻141佃々m栅ctⅡ缁or3“lprIncil∞]㈣帅i∞l ∞％orsⅢm∞mrock 加载．急居叮加大围岩损伤，窑易造成围岩疲劳破坏。 3条曲线显示．计算完毕后单元应力状忐与初始单 元班力状态相差不犬，说明地震后围岩应力状态未 产生过大扰动调整。由于地震中单元应力状态随时 闻小断变化，应力状态4；能客观反映围岩稳定情况。 笔者认为，应采用损伤程度和应力状态共同描述围 岩在地震中的劣化程度和围岩稳j；i三情况。 @从图15可以看出．随着地震渡的传播，部 分围岩长对问处于弹性应力和塑性应力空营变化状 态．在13s时，整体模型塑性损伤单元数达到峰值． 抟56545个单元。塑性损伤破坏区范田最大时捌与 地震波荷载位移峰值时刻基本一致。在I3s时刻-主 厂房边墙国岩大部进^塑性损伤破螭=区(见圈16卜 ：：} n60D0_ J1锶；剁k Ⅳl 目15整体模型m性损伤单元散时程曲拽 FIF'．15Tm-history㈣ofpI∞llc女口％ecI∞Ⅻ80f 脚。b“model 潜 目1613s时刻3‘机组段圈岩塑性损伤Ⅸ分布 FWl6Ⅸ矧bu“∞d_a“cd锄a舻。仰eatunitsec【i∞把when⋯13 局部拉裂区深大达7m，塑性损伤区深度达20m， 主厂房与尾闸室、3条M求渭隔岩破坏区贯穿，存 在失稳町能， 从附17可以看出，在地震衙载的循环加载F． 围岩损伤不断累积JJu剧，对围岩稳定造成牧人影响。 洞室用岩的塑性损伤区范围较开挖完毕时有所扩大 f见图】8)．围岩埙伤程度有显著增加．说明地震对 涧室围岩稳定肯很太影响。 目I7 i』目顶拱用柑接协系数时#∞线FiE：17■㈣histmy⋯ofdanmgeoem⋯lof⋯d1％rockattopmh㈣in 一Ⅻ10L∞ 、 - 目1825s"女性E体*室围誊破％Ⅱ*度日 Fi918Damaged叩mof“m汕曲gr