midas gts理论分析_1.doc

第一篇 MIDAS/GTS的分析功能 第一篇 MIDAS/GTS的分析功能岩土分析(geotechnical analysis)与一般的结构分析(structural analysis)有较大差异。一般的结构分析注重荷载的不确定性，所以在分析时会加载各种荷载，然后对分析结果进行各种组合，最后取各组合中最不利的结果进行设计。岩土分析注重的是施工阶段和材料的不确定性，所以决定岩土的物理状态显得格外重要。在岩土分析中应尽量使用实体单元真实模拟围岩的状态、尽量接近地模拟岩土的非线性特点以及地基应力状态(自应力和构造应力)、并且尽量真实地模拟施工阶段开挖过程，这样才会得到比较真实的结果。 优秀的岩土分析程序应能真实地模拟现场条件和施工过程，并应为用户提供更多的材料模型和边界条件，让用户在做岩土分析时有更多的选择。MIDAS/GTS不仅具有岩土分析所需的基本分析功能，并为用户提供了包含最新分析理论的强大的分析功能，是岩土和隧道分析与设计的最佳的解决方案之一。 MIDAS/GTS中提供的的分析功能如下：A. 静力分析 (static analysis)线弹性分析 (linear elastic analysis)非线性弹性分析 (nonlinear elastic analysis)弹性分析 (elastoplastic analysis)B. 施工阶段分析 (construction staged analysis)C. 渗流分析 (seepage analysis)稳定流分析 (steady state seepage analysis)非稳定流分析 (transient state seepage analysis)D. 渗流-应力耦合分析 (seepage stress analysis)E. 固结分析 (consolidation analysis)排水/非排水分析 (drained/undrained analysis)固结分析 (consolidation analysis)F. 动力分析 (dynamic analysis)特征值分析 (eigenvalue analysis)反应谱分析 (response spectrum analysis)时程分析 (time history analysis)1. 静力分析 (Static Analysis)静力分析是指结构不发生振动状态下的分析，一般来说外部荷载的频率在结构的基本周期的1/3以下时可认为是静力荷载。静力分析的类型如下：A. 线弹性分析 (linear elastic analysis)B. 非线性弹性分析 (nonlinear elastic analysis)C. 弹塑性分析 (elastoplastic analysis)1.1 线弹性分析岩土分析中的线弹性分析是将围岩材料视为线弹性，分析其在静力荷载下的响应。岩土材料的线弹性阶段仅发生在荷载加载初期应变非常小时。线弹性分析不考虑破坏将应力-应变关系理想化为直线，计算相对简单方便。从理论上说，有限元方程式的表现形式是基于虎克(Hooke)法则的线弹性方程式，非线性分析或弹塑性分析也可以按线弹性方程式的形式进行求解计算。从1990年开始，在实际设计中才开始大量使用非线性分析和弹塑性分析。其原因是非线性分析和弹塑性分析的收敛计算需要较长的时间，无论从硬件还是从软件上都还不能满足实际设计的需要。随着计算机分析速度的提高以及分析技术的发展，为非线性分析和弹塑性分析在实际设计中的应用提供了可能。但是线弹性分析以其特有的计算效率在非线性特点不是很明显的材料的分析中，作为初步分析还在大量使用。 土木领域的大部分问题可以概括为两个问题，一个是“结构在给定的荷载作用下是否安全？”，一个是“结构到完全破坏前的变形有多大？”。为了获得地基的变形需要地基的应力-应变关系，但是众所周知岩土材料的本构关系相当复杂，与材料的构成、孔隙比、应力历程以及加载方式均有关系。在实际设计中，为了便于计算会将岩土的应力-应变关系简化成一些理想化的本构关系。虽然仅用弹性模量和泊松比的变化来描述岩土特性不是很准确，但是对模拟一些特定的岩土材料还是非常有效的。在此要注意的是对弹性模量的定义。一般来说，经常使用的弹性模量包括切线模量(Tangent modulus)和割线模量(secant modulus)。完全线弹性材料的切线模量和割线模量相同，但是在岩土等非线性材料中一般使用的是所关心的应力范围内的割线模量，并将其称为变形模量(deformation modulus)。MIDAS/GTS的线弹性分析(linear static analysis)中使用的基本方程中的平衡方程式(equilibrium equation)如下。 (1.1)且: 结构物的刚度矩阵 (stiffness matrix): 位移向量 (displacement vector): 荷载向量(load vector)或不平衡力向量(unbalanced force vector)通过平衡方程式求得位移向量。这样已知荷载和刚度计算位移的方法叫位移法 (displacement method)。利用求得的位移通过变形协调方程(compatibility equation)可以得到应变，然后通过本构方程(constitutive equation)可获得应力。模型发生变形时，模型内部的任意点的坐标(x, y, z)将移动到新的坐标(x+u, y+v, z+w)位置。单元不是刚体时，位移向量(u, v, w)在单元内部是连续变化的，这种变化可以用x、y、z坐标的函数来表现。如下图所示，任意空间上分别具有微小长度dx、dy、 dz的三个具有方向的纤维(fiber)在变形后具有新的方向。图 1.1 位移(u, v, w)的定义 (1.2)在弹性材料上施加单轴应力时，将产生轴向应变。 (1.3)且 ,: x, y, z轴向应变: 弹性模量: 泊松比施加剪切应力时，剪切应变的计算公式如下。 (1.4)且，是剪切模量(shear modulus)。剪切模量与弹性模量、泊松比的关系如下。 (1.5)岩土材料的体积变形率如下： (1.6)且， (1.7)所以体积模量 K (bulk modulus) 可使用下面公式表示。 (1.8)在岩土上使用体积弹性模量K(bulk modulus)和剪切模量G(shear modulus)的概念虽然不是很准确，但是比E和n表现得更简单更明确，使用起来更方便。下图说明的是K和G的物理意义。xzconditionAccording to the loading of the stress incrementAccording to the magnitude compressionConfined compressionIsotropic Simple shearUniaxial loading=es11ddseseDDSecant modulusTangent moduluseStrainStresszszeYoungs moduluszzEse=xztxzgShear moduluszzEse=0sBulk moduluszxzxKtt=Constrained moduluszzMy图1.2 Various Types of modulus在左右边界被约束的状态下正常发生变形时，可计算出侧限模量(constrained modulus)。特别是当时，水平方向应力和侧限模量的关系如下。 (1.9) (1.10)通过现场试验可以得到上述各种弹性模量中的一个，通过适当的转换后可以应用到实际设计当中。一维固结的边界条件与计算侧限模量时的边界条件相同，所以侧限模量与软弱地基的一维固结特性密切相关。下面的表1.1中整理了侧限模量和各种一维固结特性参数的关系式。表 1.1 固结特性参数和侧限模量的关系与固结相关的参数与的关系coefficient of volume change, 体膨胀系数coefficient of compressibility, 压缩系数compression index, 压缩指数表 1.2 岩石以及其他材料的弹性模量和泊松比岩土材料弹性模量 (tonf/m2)泊松比闪岩(Amphibolite)9.412.1 1060.280.30硬石膏(Anhydrite)6.8 1060.30辉绿岩(Diabase)8.711.7 1060.270.30闪长岩(Diorite)7.510.8 1060.260.29白云石(Dolomite)11.012.1 1060.30纯橄榄岩(Dunite)14.918.3 1060.260.28含长石的片麻岩(Feldspathic gneiss)8.311.9 1060.150.20辉长岩(gabbro)8.911.7 1060.270.31花岗岩(granite)7.38.6 1060.230.27冰(ice)7.1 1060.36石灰石(limestone)8.710.8 1060.270.30大理石(marble)8.710.8 1060.270.30云母片岩(mica Schist)7.910.1 1060.150.20黑曜石(obsidian)6.58.0 1060.120.18奥长岩(oligoclasite)8.08.5 1060.29石英岩(quartzite)8.29.7 1060.120.15岩盐(rock salt)3.5 1060.25板岩(slate)7.911.2 1060.150.20铝(aluminum)5.57.6 1060.340.36钢(steel)20.0 1060.280.29表1.2中的弹性模量是采用无裂纹的小的试验体在实验室实验获得的完整岩(intact rock)的弹性模量。所以考虑现场条件，要考虑尺寸效应、岩体内的不连续性等因素应采用折减后的弹性模量。图1.3是各种岩石质量指标RQD(Rock Quality Designation)对应的弹性模量实测值图形。RQD是指10cm以上长度的岩心累计的钻孔长度比。即使RQD为100%也不能视为完整岩，但是RQD值越高，岩石品质越好。风化越严重，岩石的RQD越低。图 1.3 RQD与弹性模量折减率(EL/EM)的关系由上图可知，RQD为70%时，实验室的弹性模量就要折减20%。三维条件下，材料的应力-应变关系如下： (1.11)将上述矩阵求逆得 (1.12)且，即 (1.13) (1.14)且， 变形协调方程的矩阵如下: (1.15)且， (1.16)1.2. 非线性弹性分析岩土分析中的非线性弹性(nonlinear elastic)和弹塑性(elastoplastic)材料特性均属于材料非线性分析。所谓材料非线性是指应力与应变关系的非线性。非线性弹性材料是指材料的弹性特性随分析结果而变，其代表为像邓肯-张模型(Duncan-Chang model)这样的双曲线模型(hyperbolic model)。该模型的应力-应变关系为双曲线形状，基床系数是地基的约束 (confinement)应力和剪切应力的函数。非线性材料模型的参数可以通过三轴试验或文献中较为容易地获得，所以被应用于很多研究当中，但是其缺点是不能考虑破损后的刚度降低。图 1.4 引自: Duncan-Chang model 应力-应变曲线1.3. 弹塑性分析地基分析也可以概括为对判断在已知荷载作用下“地基具有多少安全度”的问题和“地基可以发生多大的变形”的问题。如果说线弹性分析是分析变形能力 (deformability)，则弹塑性分析则是同时分析稳定性(stability)和变形能力。地基的稳定性一般由剪切强度决定，变形能力由弹性特性和剪切特性决定。荷载作用大于地基的剪切强度时地基将产生塑性区域，随着塑性区域的发展最后达到破坏状态。但是不能说产生了塑性区域结构就一定不稳定，因为被弹性区域包围的塑性区域 (confined yield zone)不能生成破坏面，这样的局部破坏不一定会发展成为整体破坏。使用荷载作用下产生的累加位移计算得应变包括弹性应变和塑性应变。 (1.17)且，: 总应变: 弹性应变: 塑性应变在计算公式中将要使用的基本概念如下： 塑性变形的屈服标准 (yield criteria) 定义塑性变形用的流动法则 (flow rule) 变形硬化的硬化法则 (hardening rule)1.3.1. 屈服标准定义弹性区域的边界的屈服函数(或者荷载函数)F如图1.5所示。 (1.18)且，: 当前的应力: 等效(equivalent)或有效(effective)应力: 的硬化因子: 等效(equivalent)塑性应变塑性理论中屈服函数的值为正的应力状态是不存在的。产生屈服时，塑性变形逐渐累加直到屈服函数减少到零时，应力状态要不断修正。这样的过程叫塑性修正(plastic corrector)阶段或蜕化映射(return mapping)。图 1.5 关联流动准则与奇异点1.3.2 流动准则使用图1.5的流动准则定义塑性变形。 (1.19)且， : 塑性变形的方向: 定义塑性变形大小的塑性系数函数g为“塑性势能(plastic potential)”，一般使用应力不变量(stress invariant)定义。另外，塑性势能函数g与屈服函数F相同时，即g=F时称为“关联流动(associated flow)准则”，gF时称为“非关联流动(non-associated flow) 准则”。MIDAS/GTS的所有材料模型使用关联流动准则，即塑性应变向量垂直于屈服面，所以上面公式可以使用下面公式表现。 (1.20)如图1.5所示在图中角点或平面上，产生不能确定塑性流动的方向的奇异点(singular point)，对这些点需要做特殊处理。 1.3.3 本构方程标准塑性本构方程(constitutive equation)形成步骤如下。应力由应变变化率向量的弹性部分决定，即 (1.21)且， : 弹性刚度矩阵应力始终要在屈服面上，所以要满足下面的协调条件(consistency condition)。微小的应变变化率如下： (1.22)使用完全牛顿-拉普森(Newton-Raphson)迭代计算时，使用协调刚度矩阵(consistent stiffness matrix)会加快收敛速度。 (1.23)且， 1.3.4 应力积分应力积分使用显式前进欧拉方法和隐式后退欧拉方法。A. 显示前进欧拉方法(explicit forward euler algorithm with sub-incrementation) B. 隐式后退欧拉方法 (implicit backward euler algorithm) (a) 交点A的位置(b) 由A沿切向移动到C后修正到D图1.6 显式前进欧拉方法图1.7 显式前进欧拉方法的子增分(sub-incrementation)图1.8 隐式后退欧拉方法显式方法中塑性流动的方向是在交叉点，即弹性应力增量穿过屈服面的点(图1.6的A)计算的，而隐式方法是在最终应力点(图1.8的B)上计算的。显式方法相对简单且直接对应力积分，即不必在高斯点(gauss point)反复迭代计算，但也有下列缺点。 在一定条件下才能稳定。 为了满足准确度，在修正应力过程中需要子增量积分。 为了修正偏离屈服面的程度需要使用人为回归方法。另外，使用该方法不能构成协调刚度矩阵(consistent stiffness matrix)。但是隐式方法不必使用子增量或人为回归方法也可以得到较为精确的结果，并且与给定的条件无关相对稳定。但是隐式方法需要在高斯点进行反复迭代计算。使用隐式方法可以构成协调刚度矩阵，使用Newton-Raphson方法计算，也可以提高迭代计算的效率。(1) 显式前进欧拉方法的步骤Step-1 : 计算应变增量。 (1.24)且，: 应力-应变关系行列式: 位移的变化量 Step-2 : 计算假定为弹性变形的弹性应力(图1.6(a)的B点)。 (1.25)公式(1.25)和(1.26)的角标参见图1.6。Step-3 : 计算得到的应力在屈服面以内时，则完成应力修正；如果在屈服面外则根据塑性变形回归到屈服面。 Step-4 : 计算交叉应力。弹性应力的增量可分为容许应力增量和不容许的应力增量，交叉应力使用下面公式计算(参见图1.6(a)的A点)。 (1.26)Step-5 : 屈服后应力点在屈服面上移动，可使用m个不允许的应力增量近似模拟(参见图1.7)。子增量的数量与误差的大小有关，使用下面公式计算。 (1.27)Step-6 : 最终应力状态不在屈服面上时，使用人为回归方法移动到曲面上(参见图1.7的E点)。 (1.28)注意: 屈服面的形状使用各子增量的结束点使用硬化法则修正。 卸载(unloading)时假设为弹性。(2) 隐式后退欧拉法则的步骤隐式方法使用下面公式计算最终应力，角标参见图1.8。 (1.29)公式中C点是未知点，使用Newton Raphson方法反复迭代计算。任意向量为当前的应力与后退欧拉应力间的差。 (1.30)反复迭代计算的目的是将向量减少到接近于零，最终应力应满足屈标准。将假定的弹性应力按台劳(Taylor)级数展开。 (1.31)且， : 的变化量: 的变化量将上式设为零，解，得下面公式。 (1.32)将屈服函数使用台劳展开，得 (1.33)且，: 有效塑性应变。由此可得(公式1.34)，进一步可计算最终应力。 (1.34)1.3.5 非线性方程的迭代计算方法 在前面已经讲述了线性分析的有限元平衡方程式。但是当材料为非线性时，整体刚度矩阵将变成非线性，需要使用反复迭代计算方法解非线性方程式。一般来说，非线性分析就是查找荷载作用下的结构的平衡状态。在任意阶段i的平衡问题可归纳为公式(1.35)。 (1.35)且，: 节点不平衡力: 外部荷载: 由单元应力计算的内力: 节点位移迭代过程从假设的平衡状态开始分析，不平衡力()视为零，外部荷载()是已知的外部荷载，内力()从单元应力计算而得。利用公式(1.35)反复迭代计算，最后获得收敛解。解非线性方程式的方法有很多，MIDAS/GTS中提供了初始刚度法(constant stiffness method)和Newton-Raphson。(1) 初始刚度法 使用弹性初始刚度解方程式。节点位移的反复计算公式如下: (1.36)或者 (1.37)在j阶段的第i次迭代过程中计算的位移如下： (1.38)第一次生成整体刚度矩阵后计算其逆矩阵，每次迭代中仅反复计算不平衡节点力和反复使用公式(1.37)。该方法使用初始刚度矩阵，所以收敛速度较慢。但是收敛性较好且迭代速度也较快。迭代计算过程参见图1.9。图1.9 初始刚度法(Constant Stiffness Method)MIDAS/GTS中使用的是既保留了初始刚度法的收敛性好和快速的迭代计算的优点，又加快了收敛速度的由托马斯(Thomas)建议的修正的加速法(acceleration with modified Thomas method)。该方法中的第i次迭代的位移增量公式如下: (1.39)且，: 当前阶段中到第i次迭代计算总的位移增量。: 当前迭代计算的弹性位移增量。不平衡方程式如下： (1.40)且，: 第i次迭代计算时的不平衡力。: 弹塑性刚度其中项对应的是由位移得到的内力的增量。如公式(1.41)切线刚度矩阵可分为弹性和塑性成分。 (1.41)同样可将公式(1.39)分解为弹性和塑性成分。 (1.42)重新整理塑性位移增量如下： (1.43)将公式(1.43)代入(1.42)得到如下公式: (1.44)将公式(1.39)和(1.44)组合得到如下公式： (1.45)对公式(1.45)使用最小自乘法，则可得到与修正的托马斯加速因子相同的下面的公式。 (1.46)使用修正的托马斯加速法的收敛速度比一般的初始刚度法快五倍以上，但是稳定性要低一些。(2) Newton-Raphson 法Newton-Raphson法是将公式(1.35)象下面公式那样迭代计算。 (1.47)且，i是迭代次数。某j阶段的迭代过程使用的初始位移为上一阶段(j-1阶段)结束时的位移。公式(1.47)的微分表达式如下： (1.48)且，= 切线刚度矩阵。针对节点位移的迭代计算公式如下：或第j阶段的第i次迭代计算的位移结果如下：(1.49)Newton-Raphson方法中的刚度矩阵与荷载-位移曲线相切，并在每次迭代过程中均要修正。该过程参见图1.10。要注意的是该方法不能做应变软化材料的分析，应使用初始切线刚度法。为了容易收敛，MIDAS/GTS中的初始刚度始终使用弹性刚度。图1.10 Newton-Raphson MethodMIDAS/GTS中为了提高Newton-Raphson方法的收敛性和稳定性使用线搜索(line search)选项。该选项是计算最小总势能(total potential)将位移优化从而提高收敛的方法。使用该选项时首先使用Newton-Raphson 方法计算如下的位移增量。 (1.50)在此，是迭代计算中得到的切线刚度矩阵。位移修正公式如下： (1.51)且，: 修正的位移。: 前阶段结束时的位移。简单迭代方法中的常量为1。使用线搜索功能时，该常量就是变量迭代步骤长度(step length)。总势能的计算公式如下： (1.52)在此，对势能中的位移的微分与不平衡力向量相同。为了得到总势能在特定的值下最小的条件，将上面公式(1.52)重新整理如下： (1.53)为了有常解需要满足如下条件: (1.54)在此，和为常数，所以是的函数。对于初始状态，常数的计算公式如下： (1.55)在稳定状态下是正数，是负数。所以线搜索就是查找将成为零的最小的值。由于使用准确的线搜索功能的效率不是很高，一般使用将的系数取比的系数较小的值的比较宽松的线搜索。 即 (1.56)反复迭代计算，当满足上述条件时停止计算。(3) 收敛标准 MIDAS/GTS中的收敛标准有位移收敛标准、不平衡力收敛标准、不平衡能量收敛标准。迭代计算过程中满足规定的收敛标准时，就会自动进行下一步分析。三个标准均使用向量的欧几里德范数(euclidean norm)表示。向量范数是向量大小的标志，按下列公式计算。(1.57)且，: 向量的范数(norm): 向量的第i个成分: 向量中的成分数量。位移收敛标准是到第i次迭代计算中的位移增量范数与第i次迭代计算前的位移增量范数的比值作为收敛标准。(1.58)且，: 用户定义的位移收敛标准限值。: 第k次迭代计算得到的位移增量。不平衡力收敛标准是当前阶段迭代计算的不平衡力范数与当前阶段使用的外力范数的比值作为收敛标准。(1.59)且，: 用户定义的不平衡力收敛标准限值。: 第i次迭代计算得到的不平衡力范数: 当前阶段使用的外力范数2. 施工阶段分析2.1 施工阶段分析概要岩土分析一般来说是材料非线性分析，材料的非线性特性可从岩土的初始条件获得。所谓初始条件是指施工前的现场条件，也叫原场地条件。其中原场地应力最具代表性。一般来说获得原场地的应力条件后，由此可得挖掘荷载、象莫尔-库仑这样的材料的剪切强度。然后在原场地条件下按施工顺序进行全施工阶段分析。现场的实际施工阶段非常复杂也经常发生变化，施工阶段分析一般是将其简化取比较重要的施工阶段进行分析。隧道的施工阶段例子如下： 原场地应力(自重应力+构造应力)开挖第一段支护第一段 + 开挖第二段支护第二段 + 开挖第三段支护第三段 + 开挖第四段 (重复) 使用MIDAS/GTS做施工阶段分析时可以考虑的事项如下：A. 施工阶段模拟单元的添加和删除(激活和钝化)荷载的添加和删除(激活和钝化)边界条件的变化 材料特性的变化B. 地下水分析各施工阶段的稳定流分析各施工阶段的非稳定流分析D. 渗流-应力场耦合分析利用渗流分析得到的孔隙水压力进行应力分析。程序中默认单元、荷载、边界的变化均发生在施工阶段的开始步骤(first step)，所以当实际施工过程中有这些条件的变化时，要把该变化时刻定义为一个施工阶段。也就是说，结构的变化越多，要定义的施工阶段也就越多。在任意阶段添加(激活)的单元不受前面阶段作用的荷载或应力影响，也就是说新添加的单元在激活阶段时的内部应力为零。将荷载释放系数为100%的单元删除(钝化)时，钝化掉的单元的内部应力将全部分配给留下的其他单元，从而引起剩余单元的应力发生变化。与此相反，将荷载释放系数为0%的单元删除(钝化)时，钝化掉的单元的内部应力将不分配给剩余的单元。 适当调整荷载释放系数，可以调整分配给剩余单元的应力，从而可以比较真实地模拟应力释放的过程。 隧道分析中一般不是一次性完全释放挖掘单元的应力，而是随着喷锚支护阶段逐渐释放。此时可指定在不同施工阶段的荷载释放系数。MIDAS/GTS的施工阶段分析采用的是累加模型，即每个施工阶段都继承了上一个施工阶段的分析结果，并累加了本施工阶段的分析结果。也就是说上一个施工阶段中结构体系与荷载的变化会影响到后续阶段的分析结果。添加单元和荷载时，只需添加本阶段增加的单元和荷载。 2.1.1 岩土单元的添加 施工阶段分析中添加的单元的初始应力状态为零。结构上将增加相当于单元自重的荷载。添加后的岩土单元的处理方式将与结构单元相同。2.1.2 删除岩土单元 岩土单元在删除前处于受力状态。假如删除单元周边作用有荷载时，剩下的单元应通过适当的应力释放，使新生成的自由面不受应力的作用。如图1.11所示，将物体A从物体B中删除。删除前各物体的应力分别为和 ，考虑了生成该应力的所有荷载。两个物体处于平衡状态，所以为了与保持平衡荷载应由物体A作用在物体B。同样荷载应作用在物体A。因此，作用在某边界上的挖掘荷载与挖掘的单元的应力状态以及这些单元的自重相关。可定义下面公式: (1.60)且，: 应变-位移关系矩阵: 挖掘体积: 单元形函数: 岩土的容重图1.11 删除单元: 荷载分配示意图2.1.3 荷载释放系数 将三维模型简化为二维模型以及简化三维模型的施工阶段时，一般使用荷载释放系数模拟删除单元后各施工阶段的效果。MIDAS/GTS中可指定任意阶段的荷载释放系数，例如从开挖连续三个施工阶段中的应力释放比例假定为40%、30%、30%，在MIDAS/GTS中的开挖阶段中可分别定义0、1、2阶段的荷载释放系数为0.4、0.3、0.3(参见图1.12)。32图1.12 输入荷载释放系数对话框2.1.4 位移初始化(清零) 在施工阶段分析过程中有时要将位移清零。例如初始地应力作为初始的荷载条件其位移应为零。MIDAS/GTS中可以任意指定阶段做位移初始化，这样在需要事先做一些分析(比如地应力的计算、渗流)后再清零的施工阶段分析中非常实用.2.1.5 岩土材料特性的变化在施工阶段分析过程中，有时会对地基进行加固或换土处理，岩土材料有时会随时间发生硬化等。在施工阶段分析中遇到这样的问题时需要更换材料的特性. MIDAS/GTS对施工过程中修改材料特性的次数没有限制。更换材料特性对前面分析阶段的结果没有影响。352.1.6 初始自重应力初始自重应力的计算方法有两种。(1) 水平侧压力系数法给定水平侧压力系数K0按下面公式计算初始自重应力(公路隧道设计规范 JTJ D70-2004的附录J)。 (a) (b)公式中，、: 竖直方向和水平方向的初始自重地应力。 ：计算点以上第i层岩石的重度。 : 计算点以上第i层岩石的厚度。 : 计算点的孔隙水压力。在不考虑水头变化的情况下，由计算点的静水压力确定。即(为地下水 的重度，为地下水的水位差)。在计算过程中忽略剪切应力。当地面为水平面时，使用该方法没有问题。但是当地面不是水平面时，使用该方法计算得到的应力状态不能和自重形成平衡状态。为了形成平衡状态，应使用自重和上面计算的应力状态的内力差(不平衡力)进行分析。该阶段可使用没有任何变化的空阶段(null stage)进行分析。(2) 有限元法采用有限元法计算自重引起的初始应力。当地面水平时，该方法时的水平侧压力系数法的结果相同。但是当地面不是水平时，因为存在水平方向的应变，所以计算的结果与水平侧压力系数法不相同且具有剪切应力。一般来说地面不是水平面时使用该方法较好一些。但是有限元法不能计算K0大于1的情况，所以当K0大于1时可使用水平侧压力系数法计算，为了满足地应力与自重的平衡可利用空阶段(null)进行分析。3. 渗流分析(Seepage Analysis)本节中将介绍MIDAS/GTS中使用的渗流分析方法、步骤和一些技术内容。3.1 稳定流分析 稳定流分析(steady-state seepage analysis)是指岩土内部和外部的边界条件不随时间变化的分析，在分析区域内的流入量和流出量始终保持不变。非稳定流分析(transient seepage analysis)即时使用稳定的边界条件，流入量和流出量也随着时间发生变化。当存在含水土层(aquifer)，区域内边界上存在水头差(head difference)或存在流量 (flux)时，就会发生渗流现象。 渗流(seepage flow)是沿着岩土内部的粒子间孔隙流动的，一般遵循达西(Darcy)定律。在稳定流条件下通过岩土体积的渗透流量是渗透系数、水力坡降以及横断面积的之间的乘积。达西定律虽然是研究饱和状态的，但同样可适用于非饱和状态的。非饱和区域泛指从完全干燥状态到完全饱和状态的区域。饱和度在100%以下时，土壤颗粒之间的孔隙除了水以外还存在空气，饱和度非常低时水珠将以凹状附着在土壤颗粒间。 随着饱和度的降低，孔隙水压因为表面张力的影响逐渐发展成张力，所以将负的孔隙水压(negative pore pressure)也称为吸入压力(suction pressure)，一般来说饱和度越低，吸入压力越大。3.2 非稳定流分析 非稳定流分析(transient seepage analysis)的内部和外部边界随时间发生变化。 非稳定流分析与稳定流分析的区别除了边界条件随时间的变化外，非稳定流分析中使用体积含水率(volumetric water content)，以及与地下水位的变化速度密切相关的非饱和区域的含水率和孔隙率。比较土坝在干燥状态下蓄水(water filling of reservoir)与具有一定含水率状态下的蓄水时可知两种状态下土坝内部渗流达到稳定的时间有较大差异。3.3 流动法则 MIDAS/GTS中使用的流动法则为达西定律，单位面积的渗透流量公式如下： (1.61)且，: 单位面积的渗透流量: 渗透系数: 水力坡降达西定律起源于饱和土的渗透分析中，后来推广到非饱和土的渗透分析上。两种状态土的区别是非饱和土的透水系数不是常数，而是间接地随着孔隙水压的变化而变。达西法则也可以用下面公式表示。 (1.62)其中，是达西速度，水在土壤中流动时，实际平均速度是达西速度除以土的孔隙率。3.4 基本方程式在MIDAS/GTS中的流动基本方程式如下: (1.63)其中，: 总水头: x方向的渗透系数: y方向的渗透系数: z方向的渗透系数: 流量: 体积含水率: 时间该方程假定在任意位置、任意时刻微小体积的流入和流出的变化量与体积含水率的变化量相同。简单地说，x、y、z方向的流量变化与流量之和与体积含水率的变化相同。上述基本方程表现的是非稳定流的渗透方程，而稳定流状态中流入和流出量随时间没有变化，所以公式右边为零。 (1.64)体积含水率的变化依赖于应力状态的变化和土特性的变化。饱和非饱和状态的应力状态使用和两个参数表现。此处是全应力，孔隙大气压力，是孔隙水压力。MIDAS/GTS的渗流分析在假定总应力不变的情况下进行的，即土壤本身没有加载和卸载。因为没有加载，所以在非稳定流流动过程中孔隙大气压力是不变的，也就是说 是常量，对体积含水率的变化没有影响，所以体积含水率仅与的变化有关，因为是常量，所以体积含水率的变化仅与孔隙水压的变化相关。体积含水率的变化与孔隙水压变化的关系如下： (1.65)其中 是阻流系数。另外，全水头可视为压力水头和位置水头的和。 (1.66)其中，: 全水头: 孔隙水压: 水的容重: 标高重新整理公式(1.65)，得下面公式。 (1.67)将公式(1.67)代入公式(1.65)，得下面公式。 (1.68)将公式(1.68)代入公式(1.63)，得下面公式。 (1.69)标高一定时，对时间的导函数为零，变成下面的方程是。 (1.70)3.5 有限元方程使用加重残差(weighed residual)的伽辽金(Galerkin)法表现基本方程的有限元方程式如下。(1.71)其中，: 动水坡度矩阵: 单元渗透系数矩阵: 节点水头向量: 形函数向量: 单元边的单位重量: 非稳定流的阻流项: 随时间变化的水头在MIDAS/GTS的二维分析中认为所有单元的厚度都相同。所以二维有限元方程式可简化为如下公式。(1.72)其中，是单元厚度。将有限元方程用简化方式表现如下： (1.73)其中， : 单元特性军阵 = : 质量矩阵 = : 流量向量 = 公式(1.73)是非稳定流分析的基本的有限元方程式。稳定流分析的水头不随时间变化，所以MH,t项将消失，公式变成如下形式。 (1.74)3.6 时间积分非稳定流分析的有限元解是时间的函数，时间的积分可使用有限差分方法。有限元方程式使用有限差分法表达如下： (1.75)其中，: 时间增量: 0-1之间的比值。: 时间增量结束时的水头: 时间增量开始时的水头: 时间增量结束时的节点流量: 时间增量开始时的节点流量: 单元特性矩阵: 单元质量矩阵MIDAS/GTS使用后退差分法(backward difference method)，该方法的取1.0。非稳定流分析的有限元方程式的为1.0，所以 (1.75a)由公式(1.75)可知，要计算时间增量的最终阶段的水头，必须要知道开始阶段的水头。一般来说非稳定流分析必须要给出初始条件。434. 渗流-应变耦合分析(Coupled Analysis)地下水的渗透现象是由渗透区域周边边界的水头差或流量边界引起的。地下水流动过程中，因为水和土壤颗粒之间的摩擦会产生渗透力(seepage force)，由此会产生位移和应力。渗透力与对孔隙水压的积分相同。在MIDAS/GTS中利用渗流分析中计算得到的孔隙水压力(pore water pressure)计算渗透力的效果。孔隙水压是将在渗流分析中得到的全水头(total head)中减去位置水头(elevation head)而得的压力水头(pressure head)与水的容重相乘而得的。一般来说渗透力集中在全水头大小变化较大的流出边界的临近区域。这样的区域的约束压力较小，抗剪强度和抗拉强度都相对较小，所以对砂土地基做考虑渗透压的有效应力(effective stress)分析时，地基很容易破坏。所以说渗透应力耦合分析在砂土类地基的稳定性分析(analysis for stability) 中是非常重要的。4.1 有效应力地基的孔隙水压会影响总的应力。根据太沙基(Terzaghi)原理，总应力()可分为有效应力()和孔隙水压力( )。 水不能承受剪切应力，所以有效剪切应力与总的剪切应力相等。即总应力的表达式如下: (1.76)孔隙水压力可区分为稳定状态孔隙水压()和过孔隙水压()。 (1.77)稳定状态的孔隙水压可从地下水分析的结果由地下水面的高度生成。过孔隙水压由非排水材料的应力计算中产生。由此相关的过孔隙水压的计算将在后面的非排水分析中说明。4.2 耦合分析的基本方程由小位移的虎克法则可得下面的公式。 (1.78)其中是有效材料特性值。在公式(1.77)中稳定状态成分的微分为零，所以可得下面的公式。 (1.79)使用公式(1.76)和(1.79)可得下面公式。 (1.80)465. 固结分析(Consolidation Analysis)5.1 排水/非排水分析排水/非排水分析是分析渗透速度(seepage speed)和加载