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

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

乾增珍1,鲁先龙2,陈湘生3

(1.中国矿业大学力学与建筑工程学院,北京　100083;

2.国电电力建设研究所,北京　100055;3.深圳市地铁总公司,广东深圳　518000)

摘要:根据人工冻土蠕变特征和本构方程,引入经典弹塑性力学理论,将空间域和时间域进行离散化,提出了以时间增量法为基础的人工冻土有限元数值模拟计算公式应用ANSYS有限元程序分析了一个实例.并将数值模拟结果、解析解和现场实测数据进行对比分析,结果表明运用数值模拟方法可以模拟冻土的蠕变计算,揭示冻结壁变形规律,且具有一定的计算精度,从而为今后人工冻土蠕变计算提供了一种实用方法.关键词:人工冻土;蠕变;数值模拟;ANSYS程序;有限单元法中图分类号:TD265.3　　文献标识码:A

NumericalCalculationandSimulationforCreepDisplacementofArtificiallyFrozenSoil

123

QIANZeng2zhen,LUXian2long,CHENXiang2sheng

(1.SchoolofAcademics,ArchitectureandCivilEngineering,UMT,Beijing100083,China;2.ElectricPowerConstructionResearchInstitute,StatePower,Beijing100055,China;

3.ShenZhengMetroCo.Ltd,Shenzhen,Guangdong518000,China)

Abstract:Basedonthepropertiesoffrozensoilandit’sconstitutiveequationsthetimeincrement2

basedequationswereobtainedbymeansofdiscretizationoftimeandspatialdomain.ApracticalexamplewasanlaysizedbytheprogramANSYS.TheFEMsolutionsandit’scomparisonwithanalyticresultsandmeasureddataindicatethatnumericalanalysisbyprogramANSYScansimulatethecreepprocessofartificialfrozensoil,findthedeformationlawoficewallandverifythe

.Therefore,thismethodcansolvethecalculationproblemsofmathematiccreepmodeloffrozensoil

creepdisplacementinartificiallyfrozensoileffectively.Keywords:artificiallyfrozensoil;creep;numericalsimulation;programANSYS;FEM

　　自1880年德国工程师P.H.Poetsch提出人工地层冻结法(ArtificiallyGroundFreezingMethod)原理并应用于凿井工程已有120多年历史.我国首次于1955年在开滦林西风井采用冻结法凿井并获得成功.此后,冻结法在我国主要用于不稳定含水地层里的凿井工程.随着我国地下工程增多,冻结法逐渐由矿井建设工程向土木工程领域推广应用.

人工冻土最显著的特点是蠕变.在深厚表土凿井工程中,人工冻土将处于复杂的高应力状态,由冻土蠕变引起过大的蠕变变形可能引起冻结管断裂甚至引发严重的安全事故.此时蠕变变形大小决定了冻结法施工的成败,了解和掌握人工冻土蠕变特性,预计冻结壁蠕变变形大小将是设计的关键.

目前研究冻土蠕变变形的计算方法主要有理

收稿日期:20030513

作者简介:乾增珍(19752),女,安徽省宣城市人,中国矿业大学在读博士研究生,从事深基坑和边坡支护工程的优化设计以及岩土工程

的数值模拟和应用软件开发方面的研究.

　33卷27　　　　　　　　　　　　　　　　　　　　中国矿业大学学报　　　　　　　　　　　　　　　　　　第4

论计算法、数值模拟法和模型试验法.但理论计算很难全面考虑实际工程复杂的受力状况和边界条件,必须对实际问题进行简化和假设;模型试验能够在短期内获得较系统和规律性的结果,但费钱费时,而且只能获得某些特定条件下的数据;数值模拟方法具有可重复性,能较好考虑实际工程工况,已逐渐成为研究冻土蠕变特性的重要方法之一.

变,但仍处于非衰减型的蠕变变形的阶段󰂪和󰂫时,其蠕变本构方程可以表示为

=Ε

A0ΡBC

+KΡt,+()(󰃜T󰃜1)E0T

(1)

1　人工冻土蠕变特征

冻土在恒定荷载作用下发生随时间而变化的变形,即蠕变.复杂应力状态下,冻土蠕变与时间关

系曲线归纳起来,可分为两类:衰减型蠕变和非衰减型蠕变,典型的蠕变曲线如图1所示.

式中:Ρ为应力,MPa;Ε为应变,无量纲;E0(T)为冻

土单轴弹性模量,MPa;A0为试验确定的冻土蠕变常数,(MPa)-Bh-C(℃)K;B为试验确定的应力影响常数,无量纲;C为试验确定的时间影响常数,无量纲;K为试验确定的温度影响常数,无量纲;T为冻土负温值,℃;t为冻土蠕变时间,h.

在保证工程精度的前提下,考虑计算时的方便,工程中通常忽略公式(1)中第一项瞬时变形,然后将公式(1)对时间求导,得到单轴应力条件下人工冻土蠕变速率方程A0CdΕBC-1

=.KΡtdt(󰃜T󰃜+1)

(2)

图1　复杂应力状态下人工冻土典型蠕变曲线

Fig.1　Typicalcurvesofartificiallyfrozensoilunder

acombinedstressstate

图1中纵轴表示应变偏量第二不变量S2及其对时间的变化率,横轴表示时间.当应力偏量第二不变量J2较小时,蠕变呈衰减型,如图1a所示,在衰减过程中,变形速率逐渐趋近于零,蠕变变形收敛于某一个变形水平.而当J2超过某一界限值(通常称为蠕变门槛值),其蠕变呈现非衰减型,如图1b所示.非衰减型蠕变一般可分为3个阶段:非稳定蠕变阶段󰂪(蠕变速率随时间衰减),稳定蠕变阶段󰂫(蠕变速率为一个常数)和加速蠕变或破坏阶段󰂬(蠕变速率随时间增长而迅速增加,直至破坏).实际工程中,加速蠕变阶段󰂬往往经历很短时间,即可导致冻土结构物破坏.因此,对实际工程有研究意义的通常是衰减型蠕变以及非衰减型的蠕变阶段󰂪和󰂫.

但在实际工程中,受载冻结土体的蠕变计算往往处于三向应力状态,此时不能直接应用单轴应力条件下的蠕变计算公式,需要引入复杂应力状态下的蠕变本构方程.

根据经典的弹塑性力学理论[122]定义人工冻土复杂应力状态下有如下基本力学参量

1(3)[(Ρ1-Ρ2)2+(Ρ2-Ρ3)2+(Ρ3-Ρ1)2],J2=

61222

(4)[(Ε1-Ε2)+(Ε2-Ε3)+(Ε3-Ε1)],S2=

6

式中　J2为应力偏量第二不变量;S2为应变偏量第二不变量;其中Ρ1,Ρ2,Ρ3为第一、第二、第三主应力;Ε、第二、第三主应变.1,Ε2,Ε3为第一

根据研究成果[1]得到用应力偏量第二不变量J2和应变偏量第二不变量S2表示的我国人工冻土在复杂应力状态下三轴蠕变方程为

2

A0J2tS2=.

(󰃜T󰃜+1)K

B󰃗2C

(5)

定义并引入等效应力Ρi和等效应变Εi

Ρi=Εi=

233J2,3S2.

(6)(7)

将式(6)和(7)代入式(5),得到用等效应力和等效应变表达的复杂应力状态下人工冻土三轴蠕变本构方程

BC

(8)Εi=AΡit,式中　A=3-2A0(󰃜T󰃜+1)-K,其他参数同上.

对公式(8)中时间求导,得到人工冻土三轴蠕变速率公式

󰁮ΕiC-1α(9)=ACΡB.Εi=it󰁮t

B+12　人工冻土蠕变变形的数值计算

2.1　冻土蠕变的本构方程及其数值计算

根据已有研究成果[1],我国人工冻土在单轴应力条件下发生衰减型蠕变,或者虽发生非衰减型蠕

第3期　　　　　　　　　　　乾增珍等:人工冻土蠕变的数值计算及其模拟275

由上述公式可知,人工冻土蠕变过程是一个非常复杂的非线性力学过程,其应力与应变都同时随时间而变化,是时间的函数.由于蠕变变形属不可逆变形且遵循塑性理论发展,对于三维蠕变问题,通常将塑性理论推广到蠕变情况.

为了反映蠕变过程和应力应变路径有关的特性,需要得到蠕变过程中应变张量同应力张量之间的关系式.根据Prandtl2Reus塑性理论,采用张量形式表示的塑性应变增量方程如下

αp3dΕip

dΕij=Sij,

2Ρi

限个时间间隔,在每一时间间隔内假设应力水平保持不变,按前一步应力水平进行蠕变计算,这样变应力条件下求解蠕变量问题就简化成逐段常应力蠕变量的累积,即把非线性蠕变计算分段线性化,计算过程中的计算误差如图2所示.

(10)

p

式中:Εij为塑性应变张量;Sij为应力偏张量;Ρi为等

pα效应力;Ε.i为等效塑性应变

图2　数值计算中欧拉超前法计算示意图Fig.2　ThegraphforEulermethodinFEMcalculation

按时间硬化理论假定蠕变速率的各个分量和应力分量成比例有[324]

crcrαΕij=ΚSij,

i3dΕ3αcr==ΚΕij.2Ρidt2Ρi

cr从图2中可以看出,为减小误差和保持数值计算的稳定,在有限元计算时需要选用合理的时间步长,尤其是在计算的初始阶段更应如此.

(11)(12)

3　计算实例

以陈四楼煤矿主井2304m水平冻结壁蠕变位移为算例,进行冻结壁蠕变位移的数值计算.3.1　井筒概况

(11)和(12),得到人工冻土蠕变结合方程(9)、

过程中应变张量与应力张量间的关系式

3crcrB-1C-1αΕij=ACSijΡit.

2

(13)

方程(13)为有限元数值计算过程中计算蠕变应变的基本公式,式中各参数的物理意义同上.2.2　ANSYS程序中蠕变计算的实现

井筒净直径5m,冻结深度423m.据文献

[5],其蠕变关系满足蠕变方程(5),式中各参数取值分别为A0=3.31;B=3.06;C=0.39,K=1.9,根据实测平均温度T=-11℃.冻结壁内外有效半径a=3.4m,b=10.6m,掘进段高为h=2.2m.3.2　计算模型

有限元计算采用小段高(有限长厚壁圆筒)计算模型,如图3所示.

上述冻土蠕变的本构方程及其数值计算方程表明,有限元数值模拟中蠕变应变的计算虽是以经典弹塑性理论为基础,但也明显区别于一般弹塑性问题的计算,即其应力、应变都与时间密切相关.因此在有限元求解时,必须同时对时间域和空间域进行离散化,建立以时间增量法为基础的求解公式.

下面的算例采用大型通用有限元程序.该程序按欧拉超前法求解蠕应变,ANSYS计算

有如下计算算法

Εn=Εn-cr

cr

1+

crα(Ρn-1,ΕΕn-1,tn)∃t,

(14)

crcr

式中　Εn和Εn-1分别是tn和tn-1时刻的蠕变应变,Ρn-1是对应于tn-1时刻的应力水平.∃t=tn-tn-1为

计算时间步长.计算过程中不断重复上述计算过程并计算出不同时刻的位移增量,进而模拟冻土的蠕变过程,最后将计算结果加以累积就可求得蠕变变形.因此对于非线性蠕变计算,虽然蠕变过程中蠕变速率与瞬时应力水平有关,并且蠕变过程中的应力水平也逐渐变化,是时间的未知函数.但由

crα(Ρn-1,Ε于按欧拉超前法求解,求解过程中Εn-1,tn)为已知量,其计算实质是把蠕变经历的时间分成有

图3　有限段高冻结壁计算模型

Fig.3　Calculationmodeloficewall

为了研究暴露段高内的蠕变位移规律,取一个段高为研究对象,对应深度为-304～-306.2m,根据重液公式计算地压值为p=3.952MPa.模型底面取固定端约束.根据结构的对称性,取四分之一结构建立有限元模型,计算模型及其网格剖分如图4所示.

　33卷27　　　　　　　　　　　　　　　　　　　　中国矿业大学学报　　　　　　　　　　　　　　　　　　第6

图4　有限元网格剖分模型

Fig.4　TheFEMmeshedmode

图5　有限元模拟值与解析解的比较Fig.5　TheradialdisplacementsofFEMand

theoreticalsolutionsoficewall

有限元计算过程中按照用蠕变速率表示的三轴应力条件下人工冻土蠕变本构方程(9)进行,采用欧拉超前法求解公式(14)计算蠕变等效应变,最后根据式(13)求出蠕变应力和应变分量.

计算过程中蠕变计算的时间步长取0.1h.蠕变计算时间为30h.3.3　解析解

根据参考文献[5],在冻结法凿井工程中通常用下式来计算冻结壁内侧蠕变位移.

B

1--(1)(1Φ)P

A01BCB+1

,(15)U=thK1+1)(bbT1-1+a()B-1

a

a

图6　有限元模拟值与实测值的比较

Fig.6　TheradialdisplacementofFEMand

measureddataoficewall

4　结束语

人工冻土蠕变过程中应力与应变都是时间的函数,是十分复杂的非线性力学过程.有限元求解蠕变位移时,必须按时间增量法为基础建立求解公式,对时间域和空间域进行离散化,最后将计算结果加以累积求得蠕变变形.但由于对人工冻土三轴蠕变本构关系的复杂性,用有限元数值模拟方法精确地计算复杂应力状态下人工冻土三维蠕变还有待进一步深入研究.参考文献:

[1]　陈湘生.复杂应力状态下人工冻结黏土的本构关系

[A].地层冻结工程技术和应用[C].北京:煤炭工业

式中:U为冻结壁暴露段高内最大径向位移,m;a,

b为冻结壁内外有效半径,m;h为有效暴露段高,

m;P为计算处地压值,MPa,按照重液公式计算的

地压值;Φ为段高上下端约束系数,取Φ=0.4;t为空帮时间,h;其他参数同式(1),各参数代入式(15)得:U=0.011t0.393.4　计算结果与分析

有限元计算结果、按照式(15)计算的解析解和实测值[6]的相互比较分别如图5和图6所示.

从图5和6中可以看出,有限元模拟结果同解析解吻合得比较好,但是模拟数值解比实测数值大,原因可能是,有限元蠕变位移数值计算是假设一个暴露段高内土体的挖掘过程瞬间完成,然后按照蠕变求解方程(14)计算蠕变等效应变,并按照方程(13)计算蠕变应力和应变,通过累加进而计算得到暴露段高内冻结土体的蠕变位移.而实际井筒掘进过程中,开挖段高逐渐暴露,蠕变位移的发生是一个渐进的过程.因此,掘进中最大径向蠕变位移比数值模拟计算的位移小.

出版社,1995.90294.

[2]　王　仁,熊祝华,黄文彬.塑性力学基础[M].北京:

科学出版社,1998.

[3]　HarryK.

1980.

[4]　王勖成,邵　敏.

Creepanalysis[M].NewYork:Wiley,

有限单元法基本原理和数值方法

[M].北京:清华大学出版社,2001.

[5]　陈湘生.人工冻黏土力学特性研究及冻土地基离心模

型试验[D].清华大学土木水利学院,1999.

[6]　李功洲.

深井冻结壁位移实测研究[J].煤炭学报,

1995,20(1):992104.

(责任编辑　王玉浚)