（水工结构工程专业论文）地下采动对堤防危害的数值模拟.pdf

摘要 堤防是我国防洪工程体系的重要组成部分。近年来随着经济高速发展，煤 炭工业加速发展，地下采煤等人为活动对大堤等地面建筑物造成了一些破坏，影 响了堤防等正常使用。尤其是地下开采引起地表产生的位移和变形，破坏了堤防 等建筑物与地基之间的初始力学平衡状态。伴随着应力平衡的重新建立，将在堤 防等建筑物中产生附加应力，而导致大堤发生变形开裂，严重时将遭到破坏。因 此对地下采煤所造成的影响作出比较准确的预计，以预先采取足够的安全措施， 是非常重要和必要的。本文结合山东省新安二号井湖东大堤二次开采的具体情 况，采用三维快速拉格朗日有限差分法对其开采冒落进行预测，并对其上的湖东 大堤进一步进行位移及应力分析，找出影响大堤安全的危险段，危险面的位置及 具体值。为随后的工程处理提供依据。 主要内容如下： l 、回顾总结了地下采动的预测分析方法，以及我国常用的预测分析方法。 本文采用三维快速拉格朗日有限差分法进行数值模拟。 2 、简要介绍新安二号井的工程概况和已经开采工作面及其冒落后地面移动 变形的观测情况。 3 、对三维快速拉格朗日有限差分法进行了研究，并首先运用此理论对新安 二号井矿段第一阶段开采( 3 上层面) 进行预测分析，得出一系列位移及变形预 测值，并与实际观测值进行类比，验证采用n a c 3 d 的合理性。随后分析其上的湖 东大堤的应力及位移，得出第一阶段开采的危害预测。 4 、再次运用该理论对新安二号井矿段第二阶段开采( 3 下层面) 进行预测分 析，得到二次开采的综合影响结果，并综合分析二次开采对湖东大堤的应力及位 移的影响，得出影响大堤安全的危险段，危险面的位置及具体数值，完成第二阶 段开采的危害预测。 5 对本文的工作做简单的总结，并对进一步的工作提出了展望。 关键词：地下采动；损害预测分析；数值模拟；三维快速拉格朗日有限差分 法；位移应力分析 a bs t r a c t b 锄ki sa l li m p o m m tp a r ti nc l l i n e s es y s t e mo fc o n 白l l i n gj e i o o d b u tw i n lt 1 1 e l l i 曲s p e e dd e v e l o p m to fe c o n o m ya n dc o a li n d u s t n 1 0 r e 趾di n o r em a i l k i n d a c t i v i t i e s ，f o re x 锄p l e ，u n d e r g r o u n dm i n i r 培，c o l l s c i o u s l yo ru i 啪璐c i o u s l yd os o m e d 锄a g et on l eb u i l d i n 萨觚dm eb a n l ( sn e a ri t ，锄di n n u e l l c em c i rn o 锄a lu s e m o r e o v t h ed i s p l a c e m e n ta n dd e f 0 咖a t i o nb yl l i l d e 瑁的肌dm i n i n g ，h a v ed e s t r o y e di n i t i a l e q u i l i b r i 啪b e t w e e nt l l eb a r l l ( sa n dt h e 黟o u n d a d d i t i o n a ls t i c s sc 0 m e sd b o u ti nt h e b u i l d i n ga i l ds t m c t u r ew i t hm er e s u l to fr e e s t a b l i s l m e n to fe q u i l i b r i u mo nf o r c e s s y s t e i i l ，锄dm ea d d i t i o n a ls t r e s sl e a d sb u i l d i n g sa i l ds t m c t u r e st od e f l o n n ，e v e l lt ob e d e s t l d y e dw r h e l ls t a t u sb e i n gs e r i o l l s s oi ti si i t l p o r t a i l t a n de s s e i l t i a lt oa c c u r a t e l y p r e d i c tm ei i l n u e l l c eo fc o a lm i l l i n ga 1 1 d t a l ( ee n o u 曲s a f e t ym e 韶u r e si na d v a i l c e a c c o r d i n gt om es t a t i l so ft w o - s t a g em i n i n gu n d e rb a n ki nm en o 2w e l lo fx i n a n ， s h a n d o n gp r 0 v i n c c ，m i sp a p e rp r e d i c tt h ed i s p l a c e m e n t 锄dd e f 0 n n a t i o no f s u r f a c eg r o u n dw i m 凤tl a 黟a i l 舀a i la n a l y s i so fc o n t i 肌a i n3d i m e n s i o i l s ( f l a c 3 d i sa缸e e d i m e n s i o n a l e x p l i c i t f i n i t e d i f j f e r e i l c e p r o g 咖n ) ， a n d a n a l y s e m e d i s p l a c e i 】 1 e n ta 芏l ds 协鼯si i lh ud o n gb a n ka b o v et h ec o a lm i n i n g s of 如m l e 觚a l y s i sa b o v e ，w ec a nd e t e 肌i n ew h e r et h ed a n g e rp a n sa 1 1 dw h e r em ed a n g e rp l a l l e s a r e ，a n dm ea c c u r a t en u m 耐c a lv a l u e sa b o u td i s p l a c 锄e n t ，d e f o m a t i o n ，s t r e s sa i l ds o o n m a k ep r e p a r ef o rt h ef o l l o w i n ge n 垂n e e r i n gr e i n f o r c e m e n t 1 1 l cm 咖w o r kt h i sp a p e rd oa r em e s el i k eb e l o w ： 1 、t h i sp a p e rr e v i e w e da n ds u m m 撕z e dt h ep r e d i c t i o nm 武h o d sf o fu i l d e r g r o u n d m i i l i n g ，a n dm ec o m m o np r e d i c t i o nm e m o d su s e di nc h i m t 1 1 i sp a p e ri n 仃d d u c e d 缸tl 孵锄百a i la n a l y s i so fc o m i 肌血3d i m e n s i o n s ( f l a c 3 d ) m e t h o dt oi m i t a t e 2 、n l i sp a p e ri n t r d d u c e dt h eo v e i e wo fe n 舀n e 酣n gi 1 1m en o 2w e l lo f x i n a n ，a n di n t r o d u c e dt h ew o r k i n gf a c eh a v i n gb e e nm i n e da 1 1 dm e0 b s e r v a t i o no f d i s p l a c e m e n t 雠d d e f o n i l a t i o no fs u r f a c e 黟- o u n da r e rf a l l i n gb a c k w a r d 3 、h l t r c i ( 1 u c e da n ds t u d i e df l a c 3 d ，a 1 1 d1 f i r s n yu s i n gm em e o 巧t 0p r e d i c t 肌d 觚a l 弘et l l ee n d a i l g e 衄e n to fh ud o n gb a l l l ( b yt h ef i r s te x c a v a t i n gs t a g e ( 3l a y e r u p ) i nt h en o 2 、l lo fx i n a n ，锄dr e s u l t e di nas 舐e so ff o r e c a s t s s e c o n d l yv e r i 矽 觚dm a k ea i l a l o g yb e 觚e 肌m ef o r e c a s t sa i l dt h eo b s e a t i o n so fd i s p l a c e m e n ta l l d d e f 0 咖a t i o no fs u r f a c e 霉m m da r e rf 酊1 i i 培b a c 玉跚a r d a c c o r d i l l gt ot h e 删y s i so fm e d i s p l a o e m e n t 锄ds 仃e s si nh ud o n gb 锄k ，w ec a i lp r e d i c tt h e 锄d 趾g e n i l e n to f l c f i r s t s t a g eo fu i l d e 耀，o m l dm i i i i r 唱 4 、u s i n gf l a c 3 da g a i nt op 同i c ta i l d 觚a l y s e 廿l e d 觚g e n i l e n to fh ud o n g b a r l l 【b yt 1 1 es e c o r l d - s t a g ee x c a v a t i n g ( 3l a y e rd o w n ) i nm en o 2w d lo fx i n a n ， 锄dr e s u l t e di na i l o t h e rs 嘶e so ff o r e c a s t sf i n a l l y 1 1 1 e na c c o r d i n gt oa i l a l y s e so ft h e s 仃e s sa n dd i s p l a c 锄e n to ft h eh ud o n gb a l l l 【，w ec a i la l s od e t e 硼i n ew h e r et h e d a n g e rp a n s 觚dw h e r em ed a i l g e rp l 眦e sa r ei ns e c o n d - s t a g ee x c a v a t i n g ，a n dm e a c 叫a t en u m 耐c a lv a l u e s b yt h i sm e a i l s ，w ec a i la c c o m p l i s ho u r p r e d i c t i o n s 5 、s u l n m a r i z i n gt e x t u a lw o r ki nb f i e f 锄dp u t t i l l ga no u t l o o ko nm e 亿f t h e r w o l to nt b i sm e 出o d k e yw o r d s ：伽o - s t a g eu n d 哪阳u n dm i m g ；1 1 1 ep r e d i c t i o nf o re n d a n g e m l 咖； n u m 甜c a ls i i l l u l a t i o n ；f a s tl a 萨龇1 西a i l a n a l y s i so fc o n t i 肌a m3d i m e l l s i o l l s ； 趾a l y s i so fm es 懈sa l l dd i s p l a o e m 吼t 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作 及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工 作的同事对本研究所做的任何贡献均已在论文中作了明确的说明并 表示了谢意。如不实，本人负全部责任。 论文作者( 签名) ： 么轴蔓1 出 2 0 0 8 年4 月2 4 日 学位论文使用授权说明 河海大学、中国科学技术信息研究所、国家图书馆、中国学术 期刊( 光盘版) 电子杂志社有权保留本人所送交学位论文的复印件或 电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子 文档的内容和纸质论文的内容相一致。除在保密期内的保密论文外， 允许论文被查阅和借阅。论文全部或部分内容的公布( 包括刊登) 授权 河海大学研究生院办理。 论文作者( 签名) ： ! 逾五班 2 0 0 8 年4 月2 4 日 第一章绪论 1 1 开采预测分析概述 第一章绪论 1 1 1 开采后地表移动的形式【5 j 地表移动，是指采空区面积扩大到一定范围后，岩层移动发展到地表，使地 表产生移动和变形，在地表沉陷的研究中称这一过程和现象为地表移动。开采引 起的地表移动过程，受多种地质采矿因素的影响，因此，随开采深度、开采厚度、 采煤方法及煤层产状等因素的不同，地表移动和破坏的形式也不完全相同。在采 深和采厚的比值较大时，地表的移动和变形在空间和时间上是连续的、渐变的， 具有明显的规律性。当采深和采厚的比值较小( 一般小于3 0 ) 或具有较大的地质构 造时，地表的移动和变形在空间和时问上将是不连续的，移动和变形的分布没有 严格的规律性，地表可能出现较大的裂缝或塌陷坑。地表移动和破坏的形式，主 要有： ( 1 ) 地表移动盆地 在丌采影响波及到地表以后，受采动影响的地表从原有标高向下沉降，从而 在采空区上方地表形成一个比采空区面积大得多的沉陷区域。这种地表沉陷区域 称为地表移动盆地，或称下沉盆地。在地表移动盆地形成的过程中，改变了地表 原有的形态，引起了高低、坡度及水平位置的变化。因此，对位于影响范围内的 道路、建筑物、堤防等，都带来不同程度的影响。对于堤防来说，由于形成瓮地 下陷，导致堤防标高下降，从而达不到设计标准。此外，形成盆地过程中，随着 地表位移的产生，可能引起堤体的滑坡等问题，也危害了堤防的安全运行。 ( 2 ) 裂缝及台阶 在地表移动盆地的外边缘区，地表可能产生裂缝。裂缝的深度和宽度，与有 无第四纪松散层及其厚度、性质和变形值大小密切相关。若第四纪松散层为塑性 大的粘性土，一般是地表拉伸变形值超过6 l o m m m 时，地表才发生裂缝。塑 性小的沙质粘土、粘土质砂等，地表拉伸变形值达到2 3 m m m 时，地表即可 发生裂缝。地表裂缝一股平行于采空区边界发展。当采深和采厚的比值较小时， 在推进中的工作面前方地表可能发生平行于工作面的裂缝。但裂缝的宽度和深度 都比较小。这种裂缝是随工作面推进先张开而后逐渐闭合。地表裂缝的形状为楔 形，地面的开口大，随深度的增大而减小，到一定深度消失。但当采深很小、采 河海大学硕士学位论文 厚较大时，地表裂缝育可能和采空区相连通。而由地表不均匀位移引起的附近堤 体的开裂，会导致堤体的透水渗漏，特别是贯穿堤体的横向裂缝、水平裂缝以及 滑坡裂缝。它们直接影响堤体的稳定性。其中横向裂缝易发展为穿过堤身的渗流 通道，若不及时修复，可使土堤在很短的时间破坏。 1 1 2 常用的预测分析方法 1 1 2 1 基于实测资料的经验方法 基于实测资料的经验方法是通过对大量的开采沉陷实测资料的数据处理，确 定预计各种移动变形值的函数形式( 解析公式、曲线或表格) 和计算预计参数的 经验公式。在预计时，先根据丌采的地质采矿条件，运用这些经验公式求定预计 参数，再代入用上法确定的预计函数求定移动和变形值。这种方法因为是通过实 测资料来获得预计参数和函数的，所以所得出的数值比较精确。但是，这种方法 也有很多缺点：( 1 ) 实际开采沉陷只在特定的时问内发生，若不能在此时间内及 时观测到实际丌采的沉陷数据，此数据将无法在事后补测，即此丌采的实测资料 将永远无法再得到了( 2 ) 对一个尚未进行的开采可能产生的沉陷只能根据具有 类似地质采矿条件的实测资料求得的公式和参数进行预计，而测量实测资料的开 采与预计的开采的地质采矿条件常常有一定的差别，这就使没有本地区实测资料 的地区的预计的可靠性有所降低；( 3 ) 运用实测方法研究时，恶劣的天气和地理 条件将使实测变得十分困难【l j 。 1 1 2 2 影响函数法 影响函数法是介于上述经验方法和理论模拟法之间的一种方法。其实质是根 据理论研究或其他方法确定微小单元开采对岩层或地表的影响，把整个开采对岩 层或地表的影响看作采区内所有微小单元开采影响的总和，并据此计算整个开采 引起的岩层和地表的移动和变形。目前，此法所用的参数常根据实测资料求定【5 1 。 1 1 2 3 理论模拟法 理论模拟法把岩体抽象为某个数学的、力学的或数学一力学的理论模型，按 照这个模型计算受开采影响岩体产生的移动、变形和应力的分布情况。如认为岩 层和地表是一种连续的介质，则此模型属于连续介质模型；否则，就属于非连续 介质模型。该法所用的函数一般均由理论研究得出，所用的参数常用实验室试验 或理论推导求得，一般于现场实测资料没有直接关系。与丌采沉陷的实测方法相 比，模拟研究法有以下缺点：( 1 ) 岩体是一种非常复杂的介质，在模拟研究中是 不可能对它完全模拟的，这将导致模拟结果与实际结果有所差别；( 2 ) 在模拟时 2 第一章绪论 常用到一些参数，这些参数虽然可以用实验或经验的方法确定，但与岩体原始状 态的实际值常有较大的差别，这也将导致模拟结果与实际的差异。但是，它的优 点也很明显，具体有：( 1 ) 在模拟研究中，可以根据需要对某个地质条件多次进 行变换，从而研究这个地质条件与开采沉陷之间的关系；( 2 ) 在模拟研究中，可 以不受时间和地理条件的限制。 1 1 3 我国的采动预测方法 1 1 3 1 典型曲线法 典型曲线法是用无因次的典型曲线表示移动搞地主断面上的移动和变形曲 线的一种方法，它适用于矩形或近似矩形的采区的地表移动变形预计。 典型曲线法由于其分布和参数均是直接基于实测资料，所以预计误差较小， 但是典型曲线法也有它的缺点：( 1 ) 由于典型曲线法是针对某个矿区建立的，其 他矿区不能随便套用；( 2 ) 建立典型曲线要大量的实测资料；( 3 ) 典型曲线的分 布函数常常用曲线或表格给出，不便于数学处理【1 1 。 1 1 3 2 概率积分法 概率积分法是基于随机介质理论的开采沉陷预计模型。它将单元开采引起的 上覆岩层的下沉视为一随机事件，以事件发生的概率来描述岩体的沉降可能性和 沉降量，由于这种方法的基础是随机介质理论，所以又叫随机介质理论法。 概率积分法公式仅适用于矩形或近似矩形工作面的预计。该法有一重要缺 点：由于岩体本身和受采动岩体移动规律的复杂性，所采用的理论模型很难全面 地、准确地反映岩层和地表移动的规律，因此不得不进行某种假设和简化，这将 导致用其最大值和分布情况的解析公式求得的预计结果与实际情况不相符合【5 】。 1 1 3 3 负指数函数法 负指数函数法：负指数函数法是用负指数函数来表示地表下沉剖面函数的方 法，该法适用于近似矩形工作面丌采时地表移动和变形的预计。 负指数函数法比典型曲线法更便于进行数学处理和用电子计算机进行计算， 但由于其函数形式的限制，其函数值不一定能在每一点上都与实测确定的典型曲 线值相一致，求得的参数也可能与实际值有出入【5 1 。 1 1 3 4 三维快速拉格朗日法( f l a c 3 d ) 拉格朗r 单元法渊于流体力学中跟踪质团运动的一种方法，实际上是连续介 质力学中对运动的物质描述方法。它仍遵循连续介质假定，利用差分格式，按时 3 河海人学硕十学位论文 步积分求解，随构形的变化不断更新坐标，可以用于分析非线性大变形问题。 基于弹塑性理论的数值模拟方法较多，最新、最适合的方法之一是f l a c 3 d ， 它是采用三维快速拉格朗同法对连续介质进行数值分析的软件。三维快速拉格朗 日法在分析中将计算区域划分为若干个六面体单元，每个单元在给定的边界条件 下遵循指定的线性或非线性本构关系，如果单元应力使得材料屈服或产生塑性流 动，则单元网格可以随材料的变形而变形，这使它非常有利于模拟大变形问题。 它在单元剖分上采用了混合离散方法，其相邻单元的结点可以不相连接，这使得 物体的离散化更加方便；采用运动物体的动力平衡方程来模拟系统的受力变形过 程，这使得动态问题( 可以包括刚体运动) 的模拟中可以与静力学问题一样易于 解决；采用显式的求解方法使其在求解非线性本构关系时比隐式方法具有更高的 效率。 三维快速拉格朗同法的求解应用了3 种计算方法：离散模型方法将连续 介质离散为若干个互相边接的六面体单元，凡是张量( 如应力、应变) 和纯量( 如 密度) 定义在单元中，儿矢量( 如速度、位移、作用力) 定义在单元结点上； 有限差分方法变量关于时间和空间的一阶导数均用差分柬近似；动态松弛方 法应用质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。 f l a c 3 d 的主要缺点在于：求解线性问题时，其效率不如有限单元方法 ( f e m ) ；另外，其求解收敛速度取决于系统的最大固有周期和最小固有周期之 比，当单元尺寸或材料弹模相差过大时，其求解效率将降低。由于f l a c 3 d 分析 程序中引入多种理论分析模型和单元结构，因此，它可以广泛地用于岩土工程的 动力和静力稳定分析、承载力和变形计算，加固方案的设计及多体系的接触问题 的研究。又由于f l a c 3 d 是采用时间步长来进行迭代的显式有限差分程序，无论 是对动态还是静态问题，f l a c 3 d 均由显式方式求解，这使它很容易模拟动态和 大变形问题，因为对显式法来说，线性和非线的本构关系在算法上并无本质差别 f 2 8 1 。同时，由于典型曲线法只针对某一具体矿区建立的，不具备普遍性，而且以 上的预测方法，如概率积分法、负指数函数法等，只适应矿区为矩形或近似矩形， 且无法考虑具体的地质条件及施工支护条件。与之相比，f l a c 3 d 能精确的模拟 实际的地质条件和施工中采用的支护，施工的工序等，能模拟任何开采形状的矿 区，由于采用了接触单元，可以考虑矿区地质中的断裂，节理等实际情况。模拟 计算后，同时给出位移分布及应力分布等各种结果，从而可以直观的判别预测结 果。( 例如裂缝区域分布、裂缝深度等) 具有其他方法不可比拟的优越性。但是 4 第一章绪论 如果单元数量太大，或是单元剖分奇异，将使计算工作量大大增加，时间延长， 对计算机的要求较高。 1 2 问题的提出及本文研究的主要内容 随着国民经济的发展，对能源的需求也越来越大。位于苏、鲁两省交界的南 四湖( 微山湖、南阳湖、独山湖、昭阳湖) 地区是我国重要的煤炭基地之一。南 四湖地区的湖东大堤除了直接保护新安煤矿矿区等重要设施外，还承担保护4 3 8 万亩耕地，6 0 2 万人口，京沪铁路，1 0 4 国道以及微山县城的重要任务，其防洪 地位十分重要。是沂沭泗流域的防洪重点之一。但由于要穿堤入湖采煤，将使湖 东大堤产生塌陷、裂缝等破坏，直接威胁到湖东大堤的防洪安全。为了既保证华 东地区的j 下常用煤，充分利用煤炭资源，稳定矿区产量，又确保湖东大堤的防洪 安全，有必要研究预测地下采动对大堤影响的方法。 本文做了以下主要工作： 1 简要介绍了预测方法的现状及发展。 2 深入研究了地下采动的结构分析方法，三维快速拉格朗日法( f l a c 3 d ) 的基本理论。 3 介绍了新安煤矿的工程概况，和已经开采工作面及其冒落后地面沉降的观 测情况。 4 针对湖底和堤底的大面积深层二次采煤，分别运用三维快速拉格朗同法 ( f l a c 3 d ) 对地下采动，按施工顺序进行了预测模拟，首先对3 上层采煤对湖东 大堤新安矿段引起的不均匀沉降进行了计算，得出堤体及其周围地基的三维位移 场，并将结果和实测观测值进行了类比验证，进而分析了由于不均匀沉降而产生 的拉应力区和剪切破坏区域，由此评估堤体纵向裂缝、横向裂缝和剪切破坏区的 范围，从而得出了第一阶段( 3 - ：层采煤) 地下采动对大堤的影响。 5 进一步对第二阶段( 3 ，采煤) 采动引起的地面不均匀沉降进行了f l a c 3 d 的模拟和计算，并最终得出堤体及其周围地基的二次开采沉降总的三维位移场， 并分析最终在湖东大堤产生的拉应力区和剪切破坏区域，评估出大堤在的二次开 采后的纵、横向裂缝和剪切破坏区的范围，得出了最终地下采动对大堤的影响， 并对地下采动进行了一些规律性的研究。 6 提出研究中的不足及改进的可能方向。 河海人学硕十学位论文 第二章地下采动的分析方法一数值模拟( f l a c 3 d ) 2 1 三维快速拉格朗日法概述 三维快速拉格朗同法是一种基于三维显式有限差分法的数值分析方法，它可 以模拟岩石或土体及其它材料的三维力学行为。三维快速拉格朗日法在分析中将 计算区域划分为若干个六面体单元，每个单元在给定的边界条件下遵循指定的线 性或非线性本构关系，如果单元应力使得材料屈服或产生塑性流动，则单元网格 可以随材料的变形而变形，这使它非常有利于模拟大变形问题。三维快速拉格朗 日分析采用了显式有限差分格式来求解场的控制微分方程，并应用了混合单元离 散模型，可以准确模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的 弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。【7 6 】 拉格朗日单元法渊于流体力学中跟踪质团运动的一种方法，实际上是连续介 质力学中对运动的物质描述方法。它仍遵循连续介质假定，利用差分格式，按时 步积分求解，随构形的变化不断更新坐标，可以用于分析非线性大变形问题。 三维快速拉格朗同法的求解应用了3 种计算方法：离散模型方法将连续 介质离散为若干个互相边接的六面体单元，凡是张量( 如应力、应变) 和纯量( 如 密度) 定义在单元中，凡矢量( 如速度、位移、作用力) 定义在单元结点上； 有限差分方法变量关于时间和空间的一阶导数均用差分来近似；动态松弛方 法应用质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。 2 2f l a c 3 d 软件的主要特点 f l a c 3 d 是采用三维快速拉格朗同法对连续介质进行数值分析的软件。它采 用运动物体的动力平衡方程来模拟系统的受力变形过程，这使得动态问题( 可以 包括刚体运动) 的模拟中可以与静力学问题一样易于解决；在单元剖分上采用了 混合离散方法，其相邻单元的结点可以不相连接，这使得物体的离散化更加方便； 采用显式的求解方法使其在求解非线性本构关系时比隐式方法具有更高的效率。 f l a c 3 d 的主要缺点在于：求解线性问题时，其效率不如有限单元方法( f e m ) ； 另外，其求解收敛速度取决于系统的最大固有周期和最小固有周期之比，当单元 尺寸或材料弹模相差过大时，其求解效率将降低。 在f l a c 3 d 中提供了1 0 种本构模型( 表2 2 1 ) ，具有静力、动力、蠕变、渗 6 第二章地下采动的分析方法一数值模拟( f l a c 3 d ) 流和温度五种计算模式，基本涵盖了岩石、土体等弹性或非弹性材料的各种常用 本构关系，各模式间还可以相互耦合。f l a c 3 d 还提供了梁( b e 锄) 、桩( p i l e ) 、 锚( c a b l e ) 和壳( s h e l l ) 4 种结构单元( s t m c t u r a le l e m e n t ) ，用于模拟岩土工 程稳定分析中的各种支撑加固措施。为了模拟不同材料的接触面、岩体中存在的 软弱结构面等界面的作用，f l a c 3 d 中还引入了界面单元( i n t e r f a c ee l e m e n t ) 。 另外，为了便于用解决各种具体的实际问题，f l a c 如系统还提供了功能强大的 f i s h 语言和丰富的系统内部函数，以方便用户针对具体问题进行应用层面上的 二次开发。【7 6 j 正是由于f l a c 3 d 分析程序中引入上述多种理论分析模型和单元结构，因此， 它可以广泛地用于岩土工程的动力和静力稳定分析、承载力和变形计算，加固方 案的设计及多体系的接触问题的研究。 表2 2 1 模型的类型及应用范围 模型名称代表性材料戊用实例 n u l l 无材料孔洞、开挖 均一、各向同性、线性应力心变关人r t 材料( 如钢铁) 、荷载低十材 e l a s t i c 系的材料 料强度极限、计算安全系数 具有三个相互正交的弹性对称面荷载低于强度极限的柱状玄武岩 o r t h o t r o p i ce l a s t i c 的材料 t r 柚s v e r s e l yi s o t r o p i ce l 豁t i c 弹性各向异性的薄板状材料荷载低于强度极限的薄板状材料 低摩擦角材料用于吲强式有限单元程序比较的 d r u c k c r p r a g c rp l 弱t i c i 哆 模型 松散、胶结的颗粒材料；土、岩石、一般土体或岩体计算( 如边坡稳 m o h 卜c o u l o m bp l 勰t i c 毋 混凝十定、地下开挖) s t m i n - h a r d e n i n 硬化或软化时非线性的颗粒材料用于研究坡坏后的性质( 渐坍、屈 s o f t e n i n gm o h 卜c o u l o m b 服的柱体、坍方) u b i q u i t o u s - j o i n t 强度各向异性的薄板状材料层状岩层处的开挖 b 订i n e r 硬化或软化时非线性的薄板状 s t r a n h a r d e n i n s o f t e n i n g 薄板状材料的坡坏后研究 材料 u b i q u i t o u s - j o i n t 变形和剪切强度是体积改变量的粘土上的土工结构 m o d i f i e dc 锄一c l a y 函数的材料 7 河海火学硕上学位论文 2 3 三维快速拉格朗日法基本理论 显式有限差分法、运动学平衡方程、材料本构关系和连续介质的离散化是三 维快速拉格朗日法的理论基础。以下将就这几个方面对三维快速拉格朗同法的基 本理论进行阐述。为便于公式书写，公式采用了张量和标量相结合的方式表示， 并引入了爱因斯坦求和约定。【7 6 】 2 3 1 显式有限差分法 2 3 1 1 基本理论 在f l a c 3 d 中，介质中的任一点的位置、位移、速度和加速度分别用新、鳓、 和西，儡( f _ 1 ，2 ，3 ) 表示，v ，表示v f 对当前位位置坐标而求偏导。 1 物体应力应变状态的描述 给定一点的应力状态由应力张量表示，则在法向量为【n 】的面上其应力向 量i t 】可表示为如= 仃玎唧。 假定物体微元以速度【川运动，则在一个微小时段衍内，单元将产生一个微 小的应变( 称应变率孝) ，则应变率张量可表示为： 勃= ( 1 ，+ 巧，f ) 2 ( 2 - 3 一1 ) 在产生应变率张量【翻的同时，物体在该微小时段内还将产生一个微小的刚体位 移和刚体转动( 称为转动变化率) ，其转动变化率张量可表示为： = o ，- ，一吩。f ) 2 ( 2 - 3 - 2 ) 2 物体运动平衡方程 利用连续介质的冲量定理导出物体运动平衡方程 川脚警 ( 2 - 3 - 3 ) 口f 式中p 为材料密度，【6 】为单位质量的体力，平衡运动方程是物体受力单元体积 变形的控制方程。对于静力学问题，则式( 2 2 3 ) 右侧等于o ，运动平衡方程简 化为 仃+ p 抚= 0 ( 2 - 3 4 ) 3 材料本构关系 物体的运动平衡方程( 2 3 3 ) 和应变率关系等式( 2 3 1 ) 共可建立9 个等 式，而待求的未知数的个数为1 5 个，为求解还需根据材料的本构关系，建立6 个应力增量和应变增量的关系。材料的本构关系一般用下式描述： 8 第二章地下采动的分析方法一数值模拟( f l a d ) 【孑】f _ ，= h f _ ，( 盯f ，毛j ，后) ( 2 3 5 ) 其中， 孑。】是应力增量张量，旧是一个给定函数，七是一个考虑加载过程 的参数。当【民】是轮换对称的应力增量张量时，则有： 例= 譬一u 访仃万+ 仃腩u 矿 ( 2 - 3 6 ) 【伊驴j = i 矿一u 访仃矽+ 仃腩u 矽 ( 2 3 6 ) 2 3 1 2 导数的有限差分近似和介质的离散化 三维快速拉格朗同法采用了混全离散方法，连续介质区域被划分为常应变六 面体单元的集合，而在计算中又将每个六面体划分为以六面体顶点为角点的常应 变四面体单元的集合，变量均是在四面体上进行计算。六面体单元的应力、应变 值为其内四面体单元的体积加权平均。 1 用单元节点速度表示单元的应变率 如图2 3 1 所示的四面体，节点编号为l 至4 ，第n 面表示与节点n 对应的 面，设其内任一点的速率分量，则由高斯公式可得 到： j 厂矿v ，dy = j 厂sv ，刀d s ( 2 - 3 - 7 ) 其中，v 为四面体体积，s 为四面体的外表面， ，z ，为外表面单位法向向量分量。对于常应变单元，v ， 为线性分布，所在每个面上为常量，因此，由式( 2 3 7 ) 积分可得到 图2 _ 3 1 四面体单元 广壹虿u 矽w ， ( 2 - 3 - 8 ) = l 其中上标厂表示该变量为面上的变量， 速度，则有 巧= 号委v ? 1 ，= 了2 v ； 1 ，表示v ，的平均速度。对于线性变化的 ( 2 3 - 9 ) 上标，表示是节点，的变量。将式( 2 2 9 ) 代入式( 2 2 8 ) 然后两边均除以 四面体体积v ，则得到： 9 河海大学硕士学位论文 v u = 一去妻v ；刀， v u = 一万备v 彤 再将式( 2 3 1 0 ) 代入( 2 3 1 ) 则得到： 岛= 击喜( 样+ 州) ) ( 2 3 1 0 ) ( 2 3 - 1 1 ) 2 单元节点的运动平衡方程 三维拉格朗同法以节点为研究对象，将作用力和质量均集中在节点上，然后 通过运动方程在时间域上进行求解。节点的运动平衡方程，可以根据虚功原理推 导。对于静力学问题，为使系统达到平衡还需在节点运动方程中加入惯性阻尼力。 固定时间t ，则由式( 2 3 3 ) 得到该瞬时的等效静力问题平衡支配方程为： 盯玎，+ 研= 0( 2 3 1 2 ) 咖j 式中该瞬时等效体力两印仍，一j ) 。根据上述等效平衡方程，作用在四 “l 面体上节点上的力们n 应与四面体的应力和等效体力平衡，该节点力由虚功原理 推出。先对节点施加一个虚的节点速度6 【明打，则将产生一个线性虚速度场占【v 】 和虚应变率必。在该瞬时内，由内力仃，所做的功与外力( 节点力阴“和体 召】) 所做的功相等。外力所做的功可按下式求得： e = 毫啷，一+ 上以骂d y 而内力所做的功为： 卜3 ，6 乏q o y 再利用式( 2 3 11 ) 则有： i 一言( 6 v ；咿 6 v 办；。) 考虑到应力为对称张量，记z 7 = 形s d 得到： p 一 毫州t 7 将式( 2 3 1 2 ) 代入式( 2 3 1 3 ) 则外力所做的功可表示为： l o ( 2 3 1 6 ) 第一二章地下采动的分析方法一数值模拟( f l a c 3 d ) 4 e = 吖z 一+ e 6 + e 疗= l 其中矿= 以上以d y ，= 一p ，以专y 。根据前述的变分近似，在四面体 内速度场为线性变化。为便于论述，采用局部坐标系，原点位于四面体形心，坐 标轴为、x ；。虚速度可以表示为： 挑= 聊一 疗= l 3 其中的妒为线性形函数，疗= 菇+ ，打= l 以由下式得到： 疗( x ： x ： x ：) ：瓯， ( n = 1 到4 ) ，c ：为常系数，可 式中文，为n e c k e r 行列式。将局部坐标的原点定义在四面体形心后，则所有 的积分上杉d 矿均为o 。将式( 2 3 - 1 8 ) 代入式( 2 - 3 1 7 ) 则得到： 矿= 以吖懿y 疗= l 由式( 2 3 - 1 8 ) 可以解得c ；= 1 4 ，代入( 2 3 2 0 ) 则得到： 矿= 去钟以矿 同理，可以得到： 2 = 专赋沁n 删疗：l ，7 “z 将上述求得的矿和代入( 2 3 1 7 ) 则得到： e ：喜吖p + 学一，叫行警d 矿】 e = 吖弦+ 等一：，叫行等d 矿1 月= il 1 。 j 对静力学平衡问题，则有任何虚速度作用下内力虚功等于外力虚功，由式 ( 2 3 1 6 ) 和( 2 3 2 3 ) 得到 叫= 孚+ 学一上鲁刀y 协3 纠， 对四面体内的加速度场进行变分，则式( 2 3 2 4 ) 可化为： 河海人学硕士学位论文 上鲁d y = 嘲”以d y 协3 彩， 再利用局部坐标原点位于四面体形心的性质，则有： l p 一鲁d 矿= 譬( 警】一 c 2 。2 6 ， 式中丛以及前面所述的惯性阻尼均采用节点虚拟质量来代替，该虚拟质量 用来使系统的振动逐渐衰减直至平衡状态。将以上推导结果代入( 2 3 2 4 ) 就得 到： 叫= 等+ 学吲斟 协3 功， 式( 2 3 2 7 ) 即为将一个四面体单元的荷载集中到单元节点上后建立的节点 平衡方程。然后在整体坐标系下，将与某节点相关的所有单元的在该节点处的平 衡方程进行叠加，则得到节点整体平衡方程。根据牛顿定律，对每个节点则有： 矿= m 奶引d ( 乩2 3 啦) ( 2 - 3 - 2 8 ) 式中，带上标表示整体坐标下的变量，为物体节点数，只d 表示在整体坐 标下，节点，在f 方向上的不平衡力，p 表示所有与节点z 相关的单元在节点 上产生的虚拟质量之和。不平衡力f d 由下式给出： 矿= f 等+ 竿卜酽 协3 彩， 其中p 为作用在节点，相关上的外力荷载。当系统达到平衡状态时，则节点不 平衡等于0 。 3 用显式有限差分表示时间的导数 考虑到材料本构关系( 2 3 5 ) 和单元变形率与节点速度的关系( 2 3 2 9 ) ， 并考虑不平衡力的表达式( 2 3 2 8 ) ，则节点的加速可以表示如下： 鲁= 砀丢f d ( ， f 1 f 2 f 3 盱协) d 七) 卢l 到嘞。 ( 2 3 - 3 0 ) 其中 表示计算过程中所有节点的子集。在f l a c 一3 d 中，假定单元节点的速 度在一个时间段，内按线性变化，式( 2 3 3 0 ) 左边的导数用该时段的中心点处 的差分值来表示，则可得到： 霄f ( h 争= 矿( ，一争+ 嘉矿( ，pp 矿矿) 七) 1 2 第二章地下采动的分析方法一数值模拟( f l a c 3 d ) ( 2 - 3 3 1 ) 然后，节点的位置坐标利用中心差分来更新： ， ( f + ，) ：7 ( f ) + 群 o + 等) ( 2 3 3 2 ) 相应地，节点的位移按下式进行计算： ( ，+ f ) ： ( f ) + 印尸 o + 争 ( 2 3 3 3 ) 其中，取初始位移甜：( b ( o ) = o 。 通过以上推导，就可以将原计算方程中关于时间的一阶导数用代替。 4 本构模型的增量形式 在f l a c 3 d 中，假定速度在一个时间段& 内保持为常数，增量形式的本构 关系可以表示如下： 吃= 嘭( ，岛，七) ( 2 - 3 - 3 4 ) 其中，毛为轮换对称的应力增量，h ；为给定函数。 对于f 内的小位移及位移梯度可以表示为酗i f = 勺，( 勺为与结构有关的 应变改变量) 。应力增量唧= 吃+ 孑箩，其中孑孑是应力增量的修正量按下 式计算： 矛孑2 ( 盯睹缈衄一缈庸仃船) ， ( 2 3 3 5 ) 转动变化率张量的差计算公式如下： 一专善( 坩一t 拶。 ( 2 - 3 3 6 ) 以上所述提供了一种描述大变形物体的计算方法。当物体的转动张量足够 小时，则应力的修正量可以忽略不计。当物体的位移和位移梯度足够小时，则物 体的位置改变可以忽略不计，这时节点的坐标在计算中不变，这就是小变形问题。 5 节点的虚拟质量 方程( 2 3 3 1 ) 只有在物体达到平衡时才能给出合理的数值解。为使系统达 到平衡，瘵物体系统理想化为一个由线性弹簧连接的质点系统，则此质量一弹簧 系统的运动平衡方程可以用如下的矩阵形式表示。 p 宰) 一阍= 【嗍 去 其中 尸牛) 是外力， 嗣弹簧刚度，嗍为对角质量矩阵。 河海大学硕士学位论文 在采用有限差方法研究振动的质量一弹簧系统时，必须引入时间步长，这个时间 步长不能大于