人工冻土蠕变的数值计算及其模拟.pdf

第3 3 卷第3 期 2 0 0 4 年5 月 中国矿业大学学报 J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g ＆T e c h n o l o g y V o 【．3 3N o ．3 M a y2 0 0 4 文章编号1 0 0 0 1 9 6 4 2 0 0 4 0 30 2 7 3 0 4 人工冻土蠕变的数值计算及其模拟 乾增珍1 ，鲁先龙2 ，陈湘生3 1 ．中国矿业大学力学与建筑工程学院，北京1 0 0 0 8 3 ； 2 ．国电电力建设研究所，北京1 0 0 0 5 5 } 3 ．深圳市地铁总公司，广东深圳5 1 8 0 0 0 摘要根据人工冻土蠕变特征和本构方程。引入经典弹塑性力学理论，将空间域和时间域进行离 散化，提出了以时间增量法为基础的人工冻土有限元数值模拟计算公式应用A N S Y S 有限元程 序分析了一个实例．并将数值模拟结果、解析解和现场实测数据进行对比分析，结果表明运用数 值模拟方法可以模拟冻土的蠕变计算，揭示冻结壁变形规律，且具有一定的计算精度，从而为今 后人工冻土蠕变计算提供了一种实用方法． 关键词人工冻土；蠕变；数值模拟；A N S Y S 程序} 有限单元法 中围分类号T D2 6 5 ．3文献标识码A N u m e r i c a lC a l c u l a t i o na n dS i m u l a t i o nf o r C r e e pD i s p l a c e m e n to fA r t i f i c i a l l yF r o z e nS o i l Q I A NZ e n g z h e n l ，L UX i a n l o n 9 2 ．C H E NX i a n g s h e n 9 3 1 ．S c h o o lo fA c a d e m i c s ，A r c h i t e c t u r ea n dC i v i lE n g i n e e f i n g ，U M T ，B e i j i n g1 0 0 0 8 3 ，C h i n a 2 ．E l e c t r i cP o w e rC o n s t r u c t i o nR e s e a r c hI n s t i t u t e ，S t a t eP o w e r ，B e l l i n g1 0 0 0 5 5 ，C h i n a ； 3 ．S h e l lZ h e n gM e t r oC o ．L t d ，S h e n z h e n ，G u a n g d o n g5 1 8 0 0 0 ．C h i n a A b s t r a c t B a s e do nt h ep r o p e r t i e so ff r o z e ns o i la n di t ’sc o n s t i t f f t i v ee q u a t i o n st h et i m ei n c r e m e n t b a s e de q u a t i o n sw e r eo b t a i n e db ym o A l n so fd i s e r e t i z a t i o no ft i m ea n ds p a t i a ld o m a i n ．Ap r a c t m a l e x a m p l ew a sa n l a y s i z e db yt h ep r o g r a mA N S Y S ．T h eF E Ms o l u t i o n sa n di t ’sc o m p a r i s o nw i t h a n a l y t i cr e s u l t sa n dm e a s u r e dd a t ai n d i c a t et h a tn u m e r i c a la n a l y s i sb yp r o g r a mA N S Y S c a ns i m u l a t e t h ec r e e pp r o c e s so fa r t i f i c i a lf r o z e ns o i l ．f i n dt h ed e f o r m a t i o nl a wo fi c ew a l la n dv e r i f yt h e m a t h e m a t i cc r e e pm o d e lo ff r o z e ns o i l ．T h e r e f o r e ．t h i sm e t h o dc a ns o l v et h ec a l c u l a t i o np r o b l e m so f c r e e pd i s n l a c e m e n ti na r t i f i c i a l l yf r o z e ns o i le f f e c t i v e l y ． K e yw o r d s a r t i f i c i a l l yf r o z e ns o i l ；c r e e p ；n u m e r i c a ls i m u l a t i o n ；p r o g r a mA N S Y S ；F E M 自1 8 8 0 年德国工程师P ．H ．P o e t s c h 提出人工 地层冻结法 A r t i f i c i a l l yG r o u n dF r e e z i n gM e t h o d 原理并应用于凿井工程已有1 2 0 多年历史．我国首 次于1 9 5 5 年在开滦林西风井采用冻结法凿井并获 得成功．此后，冻结法在我国主要用于不稳定含水 地层里的凿井工程．随着我国地下工程增多，冻结 法逐渐由矿井建设工程向土木工程领域推广应用． 人工冻土最显著的特点是蠕变．在深厚表土凿 井工程中，人工冻土将处于复杂的高应力状态，由 冻土蠕变引起过大的蠕变变形可能引起冻结管断 裂甚至引发严重的安全事故．此时蠕变变形太小决 定了冻结法施工的成败，了解和掌握人工冻土蠕变 特性，预计冻结壁蠕变变形大小将是设计的关键． 目前研究冻土蠕变变形的计算方法主要有理 收藕日期；2 0 0 3 0 5 一1 3 作者简介；乾增珍 1 9 7 5 一 ，女，安徽省宣城市人，中国矿业大学在读博士研究生．从事深基坑和边坡支护工程的优化设计以及岩土工程 的数值模拟和应用软件开发方面的研究． 万方数据 中国矿业大学学报第3 3 卷 论计算法、数值模拟法和模型试验法．但理论计算 很难全面考虑实际工程复杂的受力状况和边界条 件，必须对实际问题进行简化和假设} 模型试验能 够在短期内获得较系统和规律性的结果，但费钱费 时，而且只能获得某些特定条件下的数据} 数值模 拟方法具有可重复性，能较好考虑实际工程工况， 已逐渐成为研究冻土蠕变特性的重要方法之一． 1人工冻土蠕变特征 冻土在恒定荷载作用下发生随时问而变化的 变形。即蠕变．复杂应力状态下，冻土蠕变与时间关 系曲线归纳起来，可分为两类衰减型蠕变和非衰 减型蠕变，典型的蠕变曲线如图1 所示． a 衰减型蠕变嘞非衰碱型蠕变 图l 复杂应力状态下人工冻土典型蠕变曲线 F i g ．1T y p i c a lC U r v e So fa r t i f i c i a l l yf r o z e ns o i lu n d e r ac o m b i n e ds t r e s ss t a t e 图1 中纵轴表示应变偏量第二不变量s 及其 对时间的变化率，横轴表示时间．当应力偏量第二 不变量山较小时，蠕变呈衰减型，如图1 a 所示，在 衰减过程中，变形速率逐渐趋近于零，蠕变变形收 敛于某一个变形水平．而当 超过某一界限值 通 常称为蠕变门槛值 ，其蠕变呈现非衰减型，如图 1 b 所示．非衰减型蠕变一般可分为3 个阶段非稳 定蠕变阶段I 蠕变速率随时f 曰衰减 ，稳定蠕变阶 段I 蠕变速率为一个常数 和加速蠕变或破坏阶 段Ⅲ 蠕变速率随时闻增长而迅速增加，直至破 坏 实际工程中，加速蠕变阶段I 往往经历很短时 间，即可导致冻土结构物破坏．因此，对实际工程有 研究意义的通常是衰减型蠕变以及非衰减型的蠕 变阶段I 和I ． 变，但仍处于非衰减型的蠕变变形的阶段I 和1 时，其蠕变本构方程可以表示为 一赤E 币7 ‰1 施。， 1 。一o r | 1 I ‘“’ 。“ 式中口为应力，M P a ；￡为应变，无量纲；E 0 为冻 土单轴弹性模量，M P a ；A 。为试验确定的冻土蠕变 常数， M P a “h _ c u c “；B 为试验确定的应力影 响常数，无量纲；c 为试验确定的时间影响常数，无 量纲；K 为试验确定的温度影响常数，无量纲} T 为冻土负温值，℃≠t 为冻土蠕变对阋，h ． 在保证工程精度的前提下，考虑计算时的方 便，工程中通常忽略公式 1 中第一项瞬时变形，然 后将公式 1 对时间求导，得到单轴应力条件下人 工冻土蠕变速率方程 塞 揣彤～． 2 出一 1 T l 十1 。o ‘ ‘ Ⅵ’ 但在实际工程中，受载冻结土体的蠕变计算往 往处于三向应力状态，此时不能直接应用单轴应力 条件下的蠕变计算公式，需要引入复杂应力状态下 的蠕变本构方程． 根据经典的弹塑性力学理论阻z 1 定义人工冻土 复杂应力状态下有如下基本力学参量 J 2 ‘} [ 巩一啦 2 口z 一如 2 d 3 一吼 2 ] ， 3 S 2 一专[ ￡1 一龟 2 e 2 ～E 3 。 e 3 一q 2 ] ， 4 式中J 为应力偏量第二不变量；S 。为应变偏量 第二不变量；其中巩，啦，口。为第一、第二、第三主应 力；E ，，￡。，岛为第一、第二、第三主应变． 根据研究成果Ⅲ得到用应力偏量第二不变量 J 和应变偏量第二不变量S 表示的我国人工冻土 在复杂应力状态下三轴蠕交方程为 z 诓 高等． ㈣ 定义并引入等效应力巩和等效应变目 厂 珥 吖3 J 2 ， 6 t 一号抵． ㈣ 将式 6 和 7 代入式 5 ，得到用等效应力和 等效应变表达的复杂应力状态下人工冻土三轴蠕 变本构方程 E ，一A 一，， 8 2 人工冻土蠕变变形的数值计算 式中A 3 一掣A 。 I TI 1 - K , 其他参数同上． 对公式 8 中时间求导，得到人工冻土三轴蠕 2 ．1冻土蠕变的本构方程及其数值计算 根据已有研究成果邮，我国人工冻土在单轴应 力条件下发生衰减型蠕变，或者虽发生非衰减型蠕 变速率公式 i 誓 A C t c ～． 9 日2 i 。 ‘‘ 9 万方数据 第3 期乾增珍等人工冻土蠕变的数值计算及其模拟 由上述公式可知，人工冻土蠕变过程是一个非 常复杂的非线性力学过程，其应力与应变都同时随 时问而变化，是时间的函数．由于蠕变变形属不可 逆变形且遵循塑性理论发展，对于三维蠕变问题， 通常将塑性理论推广到蠕变情况． 为了反映蠕变过程和应力应变路径有关的特 性，需要得到蠕变过程中应变张量闻应力张量之间 的关系式．根据P r a n d t l R e u s 塑性理论，采用张量 形式表示的塑性应变增量方程如下 2 一却 d 呓一专1 s ⋯ 1 0 ‘D 式中砖为塑性应变张量；s 。为应力偏张量；巩为等 效应力；母为等效塑性应变． 按时间硬化理论假定蠕变速率的各个分量和 应力分量成比例有0 4 1 限个时间间隔，在每一时间间隔内假设应力水平保 持不变，按前一步应力水平进行蠕变计算，这样变 应力条件下求解蠕变量问题就简化成逐段常应力 蠕变量的累积，即把非线性蠕变计算分段线性化， 计算过程中的计算误差如图2 所示． 图2 数值计算中欧拉超前法计算示意图 F i g ．2 T h eg r a p hf a rE u l e rm e t h o di nF E Mc a l c u l a t i o n 从图2 中可以看出，为减小误差和保持数值计 算的稳定，在有限元计算时需要选用合理的时间步 长，尤其是在计算的初始阶段更应如此． 豸一埚， 1 1 3 计算实例 u 3 堂3 一_r ⋯ “一2 哦出一2 毋q ‘ 、“’ 结合方程 9 、 1 1 和 1 2 ，得到人工冻土蠕变 过程中应变张量与应力张量间的关系式 弓 号A c s i 砰一1 ，一． 1 3 方程 1 3 为有限元数值计算过程中计算蠕变 应变的基本公式，式中各参数的物理意义同上． 2 2A N S Y S 程序中蠕变计算的实现 上述冻土蠕变的本构方程及其数值计算方程 表明，有限元数值模拟中蠕变应变的计算虽是U 经 典弹塑性理论为基础，但也明显区别于一般弹塑性 问题的计算，即其应力、应变都与时间密切相关．因 此在有限元求解时，必须同时对时间域和空间域进 行离散化，建立以时间增量法为基础的求解公式． 下面的算例采用大型通用有限元程序 A N S Y S 计算．该程序按欧拉超前法求解蠕应变， 有如下计算算法 ￡0 1 ， “一l ，厶一】，t 。 △f ， 1 4 式中和￡譬。分别是f 。和如一。时刻的蠕变应变， 如一，是对应于‘一，时刻的应力水平．血一k 一“一为 计算时间步长．计算过程中不断重复上述计算过程 并计算出不同时刻的位移增量，进而模拟冻土的蠕 变过程，最后将计算结果加以累积就可求得蠕变变 形．因此对于多维非线性蠕变计算，虽然蠕变过程 中蠕变速率与瞬时应力水平有关，并且蠕变过程中 的应力水平也逐渐变化，是时间的未知函数．但由 于按欧拉超前法求解，求解过程中扩 靠。，晶一，，} 。 为已知量，其计算实质是把蠕变经历的时闻分成有 以陈四楼煤矿主井一3 0 4m 水平冻结壁蠕变位 移为算例，进行冻结壁蠕变位移的数值计算． 3 ．1 井筒概况 井筒净直径5m ，冻结深度4 2 3m ．据文献 [ 5 3 ，其蠕变关系满足蠕变方程 5 ，式中各参数取 值分别为A 。 3 ．3 1 ；B 一3 ．0 6 ；c 一0 ．3 9 ，K 一1 ．9 ， 根据实测平均温度了T 一I i ℃．冻结壁内外有效 半径n 3 ．4m ，6 1 0 ．6m ，掘进段高为 2 ．2m ． 3 ．2 计算模型 有限元计算采用小段高 有限长厚壁圆筒 计 算模型，如图3 所示． 图3 有限段高冻结壁计算模型 F i g ．3 C a l c u l a t i o nm o d e lo ti c ew a l l 为了研究暴露段高内的蠕变位移规律，取一个 段高为研究对象，对应深度为～3 0 4 3 0 6 ．2m ， 根据重液公式计算地压值为户一3 ．9 5 2M P a ．模型 底面取网定端约束．根据结构的对称性，取四分之 一结构建立有限元模型，计算模型及其网格剖分如 图4 所示． 万方数据 中国矿业大学学报第3 3 卷 图4 有限元网格剖分模型 F i g ．4 T h eF E Mm e s h e dm o d e 有限元计算过程中按照用蠕变速率表示的三 轴应力条件下人工冻土蠕变本构方程 9 进行，采 用欧拉超前法求解公式 1 4 计算蠕变等效应变，最 后根据式 1 3 求出蠕变应力和应变分量． 计算过程中蠕变计算的时间步长取0 ．1h ．蠕 变计算时问为3 0h ． 3 ．3 解析解 根据参考文献E s ] ，在冻结法凿井工程中通常 用下式来计算冻结壁内侧蠕变位移． u 壶{ 杀筹} 赤甜“m s ， 式中u 为冻结壁暴露段高内最大径向位移，m } n ， b 为冻结壁内外有效半径，m ；h 为有效暴露段高， m ；P 为计算处地压值，M P a ，按照重液公式计算的 地压值；f 为段高上下端约束系数，取} 一0 ．4 ；f 为 空帮时间，h ；其他参数同式 1 ，各参数代人式 1 5 得U 一0 ．0 1 I t “” 3 ．4 计算结果与分析 有限元计算结果、按照式 1 5 计算的解析解和 实测值[ 6 2 的相互比较分别如图5 和图6 所示． 从图5 和6 中可以看出，有限元模拟结果同懈 析解吻合得比较好，但是模拟数值解比实测数值 大，原因可能是，有限元蠕变位移数值计算是假设 一个暴露段高内土体的挖掘过程瞬间完成，然后按 照蠕变求解方程 1 4 计算蠕变等效应变，并按照方 程 1 3 计算蠕变应力和应变，通过累加进而计算得 到暴露段高内冻结土体的蠕变位移．而实际井筒掘 进过程中，开挖段高逐渐暴露，蠕变位移的发生是 一个渐进的过程．因此，掘进中最大径向蠕变位移 比数值模拟计算的位移小． 5 0 4 0 目3 0 弓2 0 图5 有限元模拟值与解析解的比较 F i g ．5T h er a d i a ld i s p l a c e m e n t so fF E Ma n d t h e o r e t i c a ls o l u t i o n so fi c ew a l l f ／h 图6 有限元模拟值与实测值的比较 F i g ．6 T h er a d i a ld i s p l a c e m e n to fF E Ma n d m e a s u r e dd a t ao fi c ew a l l 4 结束语 人工冻土蠕变过程中应力与应变都是时间的 函数，是十分复杂的非线性力学过程．有限元求解 蠕变位移时，必须按时间增量法为基础建立求解公 式，对时间域和空间域进行离散化，最后将计算结 果加以累积求得蠕变变形．但由于对人工冻土三轴 蠕变本构关系的复杂性，用有限元数值模拟方法精 确地计算复杂应力状态下人工冻土三维蠕变还有 待进一步深入研究． 参考文献 E l i 陈湘生．复杂应力状态下人工冻结黏土的本构关系 L A ] ．地层冻结工程技术和应用[ c ] ．北京煤炭工业 出版社，1 9 9 5 ．9 0 9 4 ． [ 2 ] 王仁，熊棍华，黄文彬．塑性力学基础[ M ] ．北京 科学出版社，1 9 9 8 ． [ 3 ] H a r r yK ．C r e e pa n a l y s i s [ M ] ．N e wY o r k W i l e y ， 1 9 8 0 ． [ 4 ] 王勖成，邵敏．有限单元法基本原理和数值方法 [ M ] 北京t 清华大学出舨社，2 0 0 1 ． [ 5 ] 陈湘生．人工冻黏土力学特性研究及冻土地基离一I i , 模 型试验[ D ] ．清华大学土术水利学院，1 9 9 9 ． [ 6 ] 李功洲．深井冻结壁位移实测研究[ J ] ．煤炭学报， 1 9 9 5 ．2 0 】 9 9 一1 0 4 ． 责任编辑王玉浚 万方数据