《最新岩土工程数值计算ppt课件.ppt》由会员分享，可在线阅读，更多相关《最新岩土工程数值计算ppt课件.ppt（228页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、岩土工程数值计算岩土工程数值计算20112011参考文献参考文献计算土力学计算土力学 上海科学技术出版社上海科学技术出版社 朱百里朱百里土工计算机分析土工计算机分析 中国建筑工业出版社中国建筑工业出版社 龚晓南龚晓南土工数值计算土工数值计算 中国铁道出版社中国铁道出版社 钱家欢钱家欢土工原理与计算土工原理与计算 水利电力出版社水利电力出版社 钱家欢钱家欢岩土工程有限元分析理论与应用岩土工程有限元分析理论与应用 科学科学出版社出版社 谢康和谢康和Finite element analysis in geotechnical engineering David M.Potts 4 4 总刚的存储总
2、刚的存储l 一维变带宽存储一维变带宽存储由由于于整整体体刚刚度度矩矩阵阵具具有有对对称称性性、稀稀疏疏性性和和带带状状性性，采采用用变变带带宽宽下下三三角一维存储。角一维存储。DO 20 I=1,NP DO 20 J=1,JA NX=JA*(I-1)+J M=INE(I)20 MA(NX)=IWU(M,J)M=NP*JA DO 60 I=1,M K=MA(I)IF(K.LE.0)GOTO 60 DO 40 J=1,M L=MA(J)IF(L.EQ.0.OR.L.GT.K)GOTO 40 II=IDK(K)+L-K TK(II)=TK(II)+EK(I,J)40 CONTINUE 60 CONT
3、INUE将某一单元的将某一单元的KeKe集整到劲度矩阵集整到劲度矩阵KK中中l SKYLINE(轮廓线法轮廓线法)按按列列存存储储刚刚度度矩矩阵阵上上三三角角区区必必要要部部分分，按按行行存存储储下下三三角角区区中中必必要要部部分分。对对刚刚度度矩矩阵阵中中出出现现少少数数非非常常长长的的列列的的情情况况下下，存存储储要要求求不不会会剧剧烈烈增增加，很容易利用向量点积例行程序。加，很容易利用向量点积例行程序。K11 k12 k22 k13 k23 k33 k14 k24 k34DO 355 K=1,4 NRCC=3*(K-1)NR=NQ(NOD(K,IX)-1 DO 350 M=1,2+IFL
4、OW NRCC=NRCC+1 NR=NR+1 IF(ICODE.LT.2)THEN LENNC=(LOCC(NR)-LOCC(NR-1)-1)/2 DO 345 L=1,4 NCCC=3*(L-1)NCN=NQ(NOD(L,IX)-1 DO 344 N=1,2+IFLOW NCCC=NCCC+1 NCN=NCN+1 IF(NR.LT.NCN)NN=LOCC(NCN)-NCN+NR IF(NR.EQ.NCN)NN=LOCC(NR)IF(NR.GT.NCN)NN=LOCC(NR)-LENNC-NR+NCN S(NN-ISHIFT)=S(NN-ISHIFT)+C1(NRCC,NCCC)344 CON
5、TINUE 345 CONTINUE ENDIF SL(NR)=SL(NR)+ZY(NRCC)350 CONTINUE 355 CONTINUE(NCN,NR)(NR,NCN)5 5 边界条件的引入边界条件的引入l 划划 0 置置 1 法法：处理处理 ui=0 约束约束l 乘大数法：乘大数法：处理处理 ui=R 约束约束将将总总刚刚相相应应的的主主对对角角元元素素改改为为1 1，将将对对应应的的行行、列列其其它它元元素素改改为为0 0。将将荷荷载载向向量量中中相相应应的的元元素素改改为为0 0。或或预预先先将将每每个个节节点点的的方方程程编编号号，已已知知位移的节点不编号。位移的节点不编号。将
6、将总总刚刚相相应应的的主主对对角角元元素素置置一一大大数数，将将荷荷载载向向量量中中相相应应的的元元素素改改为为该大数乘该大数乘 R R。l 初等变换法：初等变换法：处理处理 两个变量有确定关系两个变量有确定关系如如轴轴对对称称问问题题，可可以以取取一一夹夹角角为为扇扇形形区区计计算算。在在斜斜对对称称ABAB边边上上的的点点只能沿只能沿ABAB上移动，所以点上移动，所以点C C上的点的位移有上的点的位移有使使用用初初等等变变换换法法消消除除一一个个相相关关方方程程，若若消消除除对对应应的的方方程程，可可以以将将所所在在的的行行乘乘tgtg 加加到到 u u 所所在在的的行行，并并将将所所在在
7、的的列列也也乘乘tgtg加加到到 u u 所在列上去。所在列上去。JP(1):NODE NUMBERJP(2):X-DIRECTION RESTRAIN INFORMATIONJP(3):Y-DIRECTION RESTRAIN INFORMATION JP(1)=I3/100 K3=I3-JP(1)*100 JP(2)=K3/10 K3=K3-JP(2)*10 JP(3)=K36 6 有限元法解题步骤有限元法解题步骤1 建立计算网格建立计算网格2 设定计算相关参数（计算精度、控制等）设定计算相关参数（计算精度、控制等）3 边界条件和内约束，力、位移、孔压等边界条件和内约束，力、位移、孔压等4
8、 选择单元类型、本构模型，并输入本构参数选择单元类型、本构模型，并输入本构参数5 计算结果提取、分析等计算结果提取、分析等岩土工程问题分析方法岩土工程问题分析方法 1.1.1 1.1 岩土工程问题控制方程的建立岩土工程问题控制方程的建立1 1 土体平衡方程土体平衡方程2 2 土体本构方程土体本构方程zxwwzpp3 3 土体几何方程土体几何方程4 4 土体有效应力原理土体有效应力原理5 5 孔隙流体平衡方程孔隙流体平衡方程6 6 渗流连续方程渗流连续方程zxw7 7 总控制方程总控制方程1.2 1.2 岩土工程基本分析方法岩土工程基本分析方法1 1 总应力分析法及其控制方程总应力分析法及其控制
9、方程总总应应力力分分析析法法与与一一般般固固体体力力学学相相同同。从从应应用用上上讲讲，一一般般用用于于不不考考虑虑渗渗流流固固结结的的情情况况，如如饱饱和和粗粗粒粒土土地地基基、透透水水性性土土料料组组成成的的土土坝坝路路堤堤的的应应力力和和变形分析以及饱和软粘土地基短期变形和稳定性分析。变形分析以及饱和软粘土地基短期变形和稳定性分析。2 2 有效应力分析法及其控制方程有效应力分析法及其控制方程在在有有效效应应力力分分析析法法中中，土土体体的的有有效效应应力力和和孔孔压压被被严严格格区区分分，并并将将土土骨骨架架变变形形与与孔孔隙隙水水的的渗渗透透同同步步考考虑虑。因因此此，有有效效应应力力
10、分分析析法法较较能能更更真真实实地地反反映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广。映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广。有有效效应应力力分分析析法法尚尚需需要要有有效效应应力力原原理理和和连连续续性性方方程程。总总应应力力法法中中只只有有位位移移变变量量且且仅仅与与空空间间有有关关，而而有有效效应应力力法法中中还还有有孔孔压压变变量量，而而且且与与空空间间和和时间均有关。时间均有关。3 3 总应力分析法总应力分析法和和有效应力分析法有效应力分析法关系关系总总应应力力法法是是有有效效应应力力法法中中当当孔孔压压p=0p=0时时的的特特殊殊形形式式。
11、在在有有效效应应力力分分析析中中如如果果采采用用与与总总应应力力分分析析相相同同的的土土工工参参数数，并并令令孔孔压压p=0p=0，所所得得结结果果即即为为总应力分析结果。总应力分析结果。总总应应力力分分析析一一般般采采用用土土体体得得不不排排水水指指标标，由由此此进进行行的的是是加加荷荷瞬瞬时时或或短短期期应应力力和和变变形形分分析析。但但也也可可采采用用土土体体的的排排水水指指标标，此此时时进进行行的的是是最最终终或或长长期期应应力力和和变变形形分分析析。当当进进行行线线弹弹性性分分析析并并采采用用排排水水指指标标，总总应应力力分分析析得得到到的结果为有效应力分析的最终（孔压消散完毕、主固
12、结完成）结果。的结果为有效应力分析的最终（孔压消散完毕、主固结完成）结果。4 4 不排水孔压计算不排水孔压计算有效应力原理有效应力原理其中其中孔压增量孔压增量不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须通过颗有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须通过颗粒接触，但接触面积很小，导致体积变化也很小，如果忽略这种体积变化。粒接触，但接触面积很小，导致体积变化也很小，如果忽略这种体积变化。则总的体积变化为则总的体积变化为对于对于饱和土体饱和土体，Ks、Kf 都比
13、土骨架模量大很多，在不考虑具都比土骨架模量大很多，在不考虑具体值时（此时准确值已不重要），可假设体值时（此时准确值已不重要），可假设Kf=Ks,得得由于由于Ks比土骨架模量大很多，如孔隙流体压缩性较大，以致比土骨架模量大很多，如孔隙流体压缩性较大，以致 Ks kf，则，则对于对于排水分析排水分析，取，取 Ke=0，加载过程中孔压不变。，加载过程中孔压不变。对于各向同性线弹性土体对于各向同性线弹性土体进行进行不排水分析时，必须设置不排水分析时，必须设置Ke。据经验，对于饱和土。据经验，对于饱和土体，只要体，只要ke足够大，土体对足够大，土体对ke的实际大小并不敏感。但的实际大小并不敏感。但Ke太
14、太大时，可能导致数值不稳定，即不排水泊松比接近大时，可能导致数值不稳定，即不排水泊松比接近0.5。专家建议专家建议设置设置 ke=Kskel，在在100和和1000之间，之间，Kskel是土骨架的体积模量。是土骨架的体积模量。u=0.1=0.3100.45200.47931000.49460.497710000.49940.49981.3 1.3 岩土工程问题的边界条件岩土工程问题的边界条件1 1 固结分析中的边界条件固结分析中的边界条件对节点位移和孔压已知的情况，可以有两种处理方法对节点位移和孔压已知的情况，可以有两种处理方法（1 1）仍给以自由度编号，在解方程时把解得的位移值用已知值置换。
15、仍给以自由度编号，在解方程时把解得的位移值用已知值置换。（2 2）已知值对应的自由度均编为已知值对应的自由度均编为0 0，把代数方程区分为：，把代数方程区分为：式中，为已知的边界值。计算中只对式中，为已知的边界值。计算中只对x x1 1建立方程，即建立方程，即Displacement+pore pressure DOFDisplacement DOF其其中中，A12x2项项在在对对各各单单元元进进行行计计算算时时即即可可求求出出。在在多多数数情情况况下下，已已知知的边界值为的边界值为0 0，上式简化为，上式简化为2 2 渗流问题的边界条件渗流问题的边界条件岩岩土土工工程程种种进进行行渗渗流流计
16、计算算的的主主要要目目的的是是求求得得渗渗流流量量和和孔孔隙隙水水头头或或水水力力坡坡降降。渗渗流流边边界界可可分分为为两两类类，即即已已知知水水头头边边界界和和已已知知流流量量边边界界，分别可以表示如下：分别可以表示如下：在渗流计算中，对于在渗流计算中，对于自由面渗流问题自由面渗流问题，其，其边界条件需特殊处理边界条件需特殊处理。(1)(1)有压渗流有压渗流有有压压渗渗流流的的计计算算域域边边界界和和边边界界条条件件是是固固定定的的，可可以以很很方方便便地地应应用用有有限限元元法法计计算算。如如图图，abab和和efef为为第第一一类类边边界界，cdcd和和fafa为为第第二二类类边边界界，
17、bcbc和和eded如如果果离离所所讨讨论论问问题题的的渗渗流流区区域域足足够够远远时时，既既可可当当作作第第一一类类边边界界，也也可可当当作作第第二二类类边边界界。设设n ni i为为计计算算域域内内部部节节点点数数，nb1nb1和和nb2nb2分分别别为为第第一一类类和和第第二二类类的的节点数，节点数，nt=ni+nb2，则渗流的有限元公式可表示为：，则渗流的有限元公式可表示为：bcdeh0h=h0h=h1h1ah=h1h=h0不透水层或或f或考虑或考虑 n=nt+nb1，并去掉下标后改写为：，并去掉下标后改写为：f fi i为为i i节节点点的的净净流流出出水水量量，对对已已知知入入渗渗
18、的的边边界界点点fifi是是给给定定的的，对对内内部部节节点点和不透水边界点均为零。和不透水边界点均为零。bcdeh0h=h0h=h1h1ah=h1h=h0不透水层或或f水头求出后，通过单元中某一截面的渗流量为：水头求出后，通过单元中某一截面的渗流量为：其中其中np为单元内节点数，为单元内节点数，lx 和和 lz为计算截面在为计算截面在x x轴和轴和z z轴的投影。轴的投影。bcdeh0h=h0h=h1h1ah=h1h=h0不透水层或或f(2)(2)自由面渗流自由面渗流在在堤堤坝坝、渠渠道道和和地地下下水水渗渗流流计计算算中中，自自由由面面或或浸浸润润面面事事先先是是未未知知的的，因此计算域本
19、身是未知的。因此计算域本身是未知的。早早期期的的研研究究工工作作是是用用试试算算法法不不断断改改变变自自由由面面位位置置，直直到到计计算算结结果果满满足足如下自由面上的边界条件为止。如下自由面上的边界条件为止。(1)(2)然然后后用用新新坐坐标标进进行行下下一一轮轮计计算算，直直到到自自由由面面上上所所有有节节点点满满足足 h-z h-z pfi 流量边界Pfbpfi 孔压边界盾构隧道施工过程分析盾构隧道施工过程分析 9.E D C B APost shield lossshield lossface lossDirection of driveA cutting head and hoodB
20、 shieldC tailD unground liningE grouted lining一、隧道施工的地面响应一、隧道施工的地面响应地层损失：地层损失：由于地层向开挖面移动，额外被开挖的超过隧道体积的土体，通常用隧道体积的百分比表示。1.Gap methodgapDmDl二、隧道开挖的模拟二、隧道开挖的模拟Rowe etal 1983Initial tunnelpositionFinal tunnelpositionThe gap=Dm-D1 at the crown,and is always fixed to zero at the invert.在有限元网格中预先预先定义空隙，代表预
21、期的地层损失。通过改变空隙的尺寸，可以反映不同的隧道施工方式。空隙在隧道周围设置，因而确定了开挖前的土体边界。将隧道仰拱设在下卧土体上，在隧道拱顶规定空隙参数。gapDmDlInitial tunnelpositionFinal tunnelpositionThe gap=Dm-D1 at the crown,and is always fixed to zero at the invert.分析时，去掉开挖面周边的剪力，并监测节点位移。当节点位移达到空分析时，去掉开挖面周边的剪力，并监测节点位移。当节点位移达到空隙合拢、土体与预先定义的衬砌位置接触时，在这些节点激活土体隙合拢、土体与预先定义
22、的衬砌位置接触时，在这些节点激活土体衬砌相衬砌相互作用。互作用。2.Convergence-confinement methodinitialLining installation/unloadingExcavation/unloadingPanet and Guenot 1982设定衬砌施工前的卸载比例，因此体积损失是一个预测值。在隧道边界节点施加力矢量（1-）F0（F0 相当于土体初始应力0）初始等于0，以后逐渐增加到1 以模拟开挖过程。在规定值d时，安装衬砌，这时边界应力减小d 0，剩余应力部分新增了衬砌应力。Ee3.Progressive softening methodinitial
23、Excavation,Lining installationand unloadingStiffness reductionand unloadingE0E0Swoboda 1979通过给土体刚度乘一降低因子降低因子，使掘进头里的土体软化。当开挖力作用于未来隧道边界时，土体出现明显软化。当隧道采用台阶法施工时，每级施工都可用上述方法。4.Volume loss methodVsVs=volume of settlement troughzVL=Ve-VtVe=excavated tunnel volumeVt=final tunnel volume类似于收敛限制法，但不规定衬砌施工前的卸载比例
24、，而是规定开挖完成后的体积损失。此法用于体积损失可以预期的给定的掘进法Forces imposedBy soil to be excavatedForces imposed ontunnel boundary due to excavation F-F 计算由被开挖土体施加于隧道边界的等效节点力F0，它被线性分配到开挖发生过程的n个增量步，即F=F0/n。在开挖的每一增量步在开挖边界施加等值反向的力矢量-F。监测每一步的体积损失，在一预想的体积损失增量步施工衬砌。衬砌施工后，对于剩下的增量步，仍然在开挖周边施加载荷-F，因而在衬砌内引起初始应力。后续过程还有可能发生进一步的地层损失，这取决于衬
25、砌的刚度。因此，在小于预期地层损失的增量步安装衬砌，完成开挖后可以得到预期的地层损失。Forces imposedBy soil to be excavatedForces imposed ontunnel boundary due to excavation F-F开挖面释放荷载的计算开挖面释放荷载的计算在盾构隧道施工过程中，由于土压舱内土体压力不能完全平衡开在盾构隧道施工过程中，由于土压舱内土体压力不能完全平衡开挖面土压力，使得土体开挖面前方土体发生趋近或者远离盾构开挖面挖面土压力，使得土体开挖面前方土体发生趋近或者远离盾构开挖面的三维位移。的三维位移。现有的模拟方法是每向前开挖一步，通过
26、开挖面处的土体原位土现有的模拟方法是每向前开挖一步，通过开挖面处的土体原位土压力与土压舱支护压力的差值压力与土压舱支护压力的差值 P 计算作用在开挖面结点上的等效节计算作用在开挖面结点上的等效节点力。这种方法的开挖面前方土体主要发生点力。这种方法的开挖面前方土体主要发生水平方向的位移水平方向的位移,而实际上，而实际上，由于开挖面卸荷引起的前方由于开挖面卸荷引起的前方土体移动是三维土体移动是三维的的，因此，该处理方法并，因此，该处理方法并不合理。不合理。本文采用基于本文采用基于 Mana 法法 的开挖释放荷载计算方法如图所示。将开的开挖释放荷载计算方法如图所示。将开挖面正前方一层土体单元作为开挖
27、卸荷单元，用该层开挖卸荷单元的收挖面正前方一层土体单元作为开挖卸荷单元，用该层开挖卸荷单元的收缩变形反映开挖引起的土体位移。缩变形反映开挖引起的土体位移。Mana 法的表达式为法的表达式为由于盾构施工中开挖面上应力由于盾构施工中开挖面上应力没有完全释放没有完全释放，因而不能像模拟基坑，因而不能像模拟基坑开挖那样，用上一步得到的单元应力直接计算新开挖面处的开挖释放荷开挖那样，用上一步得到的单元应力直接计算新开挖面处的开挖释放荷载，而应该用单元载，而应该用单元应力改变量应力改变量计算开挖释放荷载。计算开挖释放荷载。本文认为由于刀盘的搅动开挖面前方土体受到很大的扰动，可近似本文认为由于刀盘的搅动开挖
28、面前方土体受到很大的扰动，可近似认为土体由原来的的应力状态变为各向等压状态，并且压力等于土压舱认为土体由原来的的应力状态变为各向等压状态，并且压力等于土压舱支护压力支护压力 ps，土体单元应力改变量，土体单元应力改变量同时将开挖卸荷单元的模量降低。同时将开挖卸荷单元的模量降低。x，y，z=ps-x，ps-y，ps-z分别运用直接由分别运用直接由 P 计算等效节点力的方法和本文推荐方法，计算得计算等效节点力的方法和本文推荐方法，计算得到纵断面上由开挖面卸荷引起的盾构周围土体位移场如图到纵断面上由开挖面卸荷引起的盾构周围土体位移场如图 2 所示（计算所示（计算中采用的支护压力均小于开挖前土体初始应
29、力），可以看出本文推荐的中采用的支护压力均小于开挖前土体初始应力），可以看出本文推荐的方法能够较好的模拟开挖面前方土体的三维移动。方法能够较好的模拟开挖面前方土体的三维移动。刀盘超挖间隙和盾尾空隙的模拟刀盘超挖间隙和盾尾空隙的模拟为为了了使使盾盾构构能能够够顺顺利利前前行行而而减减小小盾盾构构壳壳上上的的摩摩擦擦阻阻力力，通通常常盾盾构构的的刀刀盘盘外外径径要要大大于于盾盾构构壳壳的的外外径径，从从而而在在盾盾构构壳壳外外围围产产生生一一定定厚厚度度的的间间隙隙 ，盾盾构构机机由由于于自自重重发发生生下下沉沉到到底底部部，在在横横断断面面上上形形成成超超挖挖间间隙隙，随随后后，周周 围围 土土
30、 体体 由由 于于 填填 充充 超超 挖挖 间间 隙隙 而而 产产 生生 指指 向向 盾盾 构构 内内 部部 的的 径径 向向 位位 移移当盾构前行时，由于不能做到完全同步注浆，而在盾尾产生建筑空当盾构前行时，由于不能做到完全同步注浆，而在盾尾产生建筑空隙，如图隙，如图 4 所示。图中所示。图中 t 为盾构壳的厚度，空隙所对应的衬砌管片由于为盾构壳的厚度，空隙所对应的衬砌管片由于受到前后衬砌螺栓的固定作用不能下沉，因而形成图示的盾尾建筑空隙。受到前后衬砌螺栓的固定作用不能下沉，因而形成图示的盾尾建筑空隙。模拟上述刀盘超挖间隙和盾尾建筑空隙时，首先将相应位置的盾模拟上述刀盘超挖间隙和盾尾建筑空隙
31、时，首先将相应位置的盾构壳单元和衬砌单元用刚度极小的软材料代替（活化减退），给外围土构壳单元和衬砌单元用刚度极小的软材料代替（活化减退），给外围土体施加如图所示的已知径向位移，当土体接触到盾构壳或者衬砌外壁时，体施加如图所示的已知径向位移，当土体接触到盾构壳或者衬砌外壁时，再将盾构壳和衬砌单元重新激活，使它们与土体共同作用再将盾构壳和衬砌单元重新激活，使它们与土体共同作用盾尾注浆的模拟采用对盾构壳外围的土体单元施加远离隧道中心盾尾注浆的模拟采用对盾构壳外围的土体单元施加远离隧道中心的结点荷载，其数值由注浆压力计算得到的结点荷载，其数值由注浆压力计算得到盾尾注浆的模拟盾尾注浆的模拟 F盾构法隧道
32、施工的精细模拟，张海波，岩土力学刚度迁移法是将盾构前行看成刚度和载荷的迁移过程。图刚度迁移法是将盾构前行看成刚度和载荷的迁移过程。图2是盾是盾构前行时刚度迁移法示意图。在盾首和盾尾均设有预设单元，开挖面构前行时刚度迁移法示意图。在盾首和盾尾均设有预设单元，开挖面推进时，盾首逐渐深入，预设单元刚度逐渐增加推进时，盾首逐渐深入，预设单元刚度逐渐增加;盾尾逐渐脱出，预盾尾逐渐脱出，预设单元的刚度逐渐减少。当推进了一段行程后，盾首预设单元变为外设单元的刚度逐渐减少。当推进了一段行程后，盾首预设单元变为外壳刚度，盾尾的变为空气单元，盾构前行了一个行程。壳刚度，盾尾的变为空气单元，盾构前行了一个行程。刚度
33、迁移法原理刚度迁移法原理盾构前行的同时，盾构附属的其他结构包括载荷也必须前行。盾构前行的同时，盾构附属的其他结构包括载荷也必须前行。加强环迁移与预设单元前行相同，载荷的迁移按加强环的刚度迁移加强环迁移与预设单元前行相同，载荷的迁移按加强环的刚度迁移来分配。当前加强环由前位置迁移到后位置时，载荷也随着迁移。来分配。当前加强环由前位置迁移到后位置时，载荷也随着迁移。最后，前盾构前行时，加强环、其上安装的设备及其载荷加强环成最后，前盾构前行时，加强环、其上安装的设备及其载荷加强环成为空气单元时，其上的载荷也为零。为空气单元时，其上的载荷也为零。混合刚度接触单元混合刚度接触单元盾构推进过程的有限元模拟
34、必须解决以下问题盾构推进过程的有限元模拟必须解决以下问题:(l)盾构运动。盾构在推进力作用下向前运动，但是单元网格是业已划盾构运动。盾构在推进力作用下向前运动，但是单元网格是业已划好的，仅用盾构与土体之间的接触单元错动来模拟盾构前行，与盾构好的，仅用盾构与土体之间的接触单元错动来模拟盾构前行，与盾构前行不符，见图前行不符，见图1。刚度迁移法(2)运动接触点对。盾构前行产生很大的相对位移，盾构对运动接触点对。盾构前行产生很大的相对位移，盾构对B点产生支承，点产生支承，支承面支承面BB，AB成为自由面，失去盾构的支承，成为自由面，失去盾构的支承，A不再与不再与A点形成结点对。点形成结点对。(3)力
35、学性质改变。盾构与土体间的滑动过程是有摩擦滑动，滑动在力学性质改变。盾构与土体间的滑动过程是有摩擦滑动，滑动在AB间产生，但当间产生，但当A了已经脱出盾尾时，此间不会产生摩擦。了已经脱出盾尾时，此间不会产生摩擦。随着盾构迁移，盾首和盾尾与土体的接触点在移动。以盾尾为例，随着盾构迁移，盾首和盾尾与土体的接触点在移动。以盾尾为例，接触点对接触点对AB运动的结果是最后一层，盾尾单元一部分是外壳刚度，运动的结果是最后一层，盾尾单元一部分是外壳刚度，另一部分则是空气单元刚度，两种不同的材料填充在一个单元内，见另一部分则是空气单元刚度，两种不同的材料填充在一个单元内，见图。这样的特殊单元是不能直接用图。这
36、样的特殊单元是不能直接用Goodman单元来模拟的，必须加以单元来模拟的，必须加以修正。修正。盾尾两种介质的混合单元在盾构迁移过程中，盾首和盾尾单元及其接触单元内一部分由土在盾构迁移过程中，盾首和盾尾单元及其接触单元内一部分由土体材料变为外壳材料，或由外壳材料变为空气材料，见图体材料变为外壳材料，或由外壳材料变为空气材料，见图3。单元内两。单元内两种材料影响到应力一应变矩阵，即单元刚度。假设己知单元结点的刚种材料影响到应力一应变矩阵，即单元刚度。假设己知单元结点的刚度，那么单元内任意一点刚度可以用加权形函数插值求得。度，那么单元内任意一点刚度可以用加权形函数插值求得。(l)空间空间Goodma
37、n单元内两种材料下任一点权刚度单元内两种材料下任一点权刚度式中：式中：i=Ai/A为为Ei面积权系数，反映盾尾位置对单元的影响面积权系数，反映盾尾位置对单元的影响(2)盾构外壳单元内两种材料的任一点权弹性常数盾构外壳单元内两种材料的任一点权弹性常数式中：式中：i=Ai/A为为Ei面积权系数，反映盾尾位置对单元的影响面积权系数，反映盾尾位置对单元的影响从盾构法隧道施工过程的整体来看，实际上是一种从盾构法隧道施工过程的整体来看，实际上是一种稳态掘进稳态掘进的过程，的过程，在每一时刻，盾构系统与周围土层均处于在每一时刻，盾构系统与周围土层均处于平衡状态平衡状态。然而，从。然而，从局部局部来看，来看，
38、特别是盾构切削土体这一过程，是一种复杂地特别是盾构切削土体这一过程，是一种复杂地动态相互作用动态相互作用过程。因此，过程。因此，要研究盾构机的掘进性能，必然要考虑到盾构掘进过程中的宏观平衡特要研究盾构机的掘进性能，必然要考虑到盾构掘进过程中的宏观平衡特性以及局部动态特性。性以及局部动态特性。目前，国内外对盾构目前，国内外对盾构刀盘切削土体刀盘切削土体过程的数值模拟研究非常少，过程的数值模拟研究非常少，这主要由于：一方面，要实现刀盘对土体的切削模拟，必须能够解决刀这主要由于：一方面，要实现刀盘对土体的切削模拟，必须能够解决刀盘与土体的盘与土体的接触行为接触行为模拟，处理由土体大变形、破坏等带来的
39、数值稳定模拟，处理由土体大变形、破坏等带来的数值稳定性；另一方面，刀盘对土体切削过程，实际上是一个性；另一方面，刀盘对土体切削过程，实际上是一个 3 尺度问题尺度问题，即，即，刀齿切削土体过程为厘米级别，刀盘与土体相互作用过程为米级别，整刀齿切削土体过程为厘米级别，刀盘与土体相互作用过程为米级别，整个地层环境为百米级别以上。因此，要对这一个多尺度问题进行三维动个地层环境为百米级别以上。因此，要对这一个多尺度问题进行三维动态模拟，既需要合适的算法支持，更需要强大的计算能力。态模拟，既需要合适的算法支持，更需要强大的计算能力。盾构掘进离散方法沈建奇，盾构掘进过程数值模拟方法研，上海交通大学博士学位
40、论文，2009泥水盾构掘进过程的逐步离散，是将盾构掘进的泥水盾构掘进过程的逐步离散，是将盾构掘进的连续连续稳态过程稳态过程离散离散成以一环衬砌距离为成以一环衬砌距离为掘进步长掘进步长的的逐步逐步稳态掘进过程稳态掘进过程。将泥水盾构每掘进。将泥水盾构每掘进一环衬砌距离，即认为盾构一环衬砌距离，即认为盾构-土层土层-构筑物系统处于构筑物系统处于外部和内部平衡状态外部和内部平衡状态，将这一平衡状态作为一载荷步进行迭代求解。将这一平衡状态作为一载荷步进行迭代求解。在盾构掘进过程中，相应的土体单元被在盾构掘进过程中，相应的土体单元被“杀死杀死”，同时相应位置的衬，同时相应位置的衬砌、注浆、机身单元被砌、
41、注浆、机身单元被“激活激活”，相应的荷载条件、边界条件被施加，并，相应的荷载条件、边界条件被施加，并在每一掘进步内，计算盾构在每一掘进步内，计算盾构-土体土体-构筑物系统的应力、应变，最后得到构筑物系统的应力、应变，最后得到以掘进距离为横坐标的位移、变形、应力、姿态等重要数据变化曲线。以掘进距离为横坐标的位移、变形、应力、姿态等重要数据变化曲线。下图表达了盾构机头进入土层，然后连带拖车逐渐掘进，并在掘进的过程下图表达了盾构机头进入土层，然后连带拖车逐渐掘进，并在掘进的过程中完成衬砌的拼装及安装，最后进行注浆等较为完整的盾构施工次序。中完成衬砌的拼装及安装，最后进行注浆等较为完整的盾构施工次序。
42、图图(a)为盾构开挖第一环距离示意图，此时相应位置为盾构开挖第一环距离示意图，此时相应位置土体单元土体单元被被“杀杀死死”以模拟土体被挖掘过程，同时以模拟土体被挖掘过程，同时“激活激活”预设预设机身壳单元机身壳单元，机身与土体，机身与土体接接触单元触单元随之被随之被“激活激活”。2.3.7 机身坡度及刀盘超挖模拟机身坡度及刀盘超挖模拟实际开挖过程中，盾构机头直径大于机尾直径，机身普遍具有一定坡实际开挖过程中，盾构机头直径大于机尾直径，机身普遍具有一定坡度。通过在预设度。通过在预设机身机身壳壳单元上单元上施加位移施加位移来模拟该坡度，并且以自适应来模拟该坡度，并且以自适应接触单接触单元元建立机身
43、与土体的相互作用，近似模拟实际开挖过程中，由于机身坡度的建立机身与土体的相互作用，近似模拟实际开挖过程中，由于机身坡度的存在，机身周围土体不断向隧道轴线方向偏移的复杂力学过程。存在，机身周围土体不断向隧道轴线方向偏移的复杂力学过程。因此，在因此，在“激活激活”机身单元的同时，在相应位置机身壳单元环向节点上机身单元的同时，在相应位置机身壳单元环向节点上施加径向位移施加径向位移dr，综合模拟盾构超挖及机身坡度，综合模拟盾构超挖及机身坡度.图图(c)为盾构完全进入土层后，超挖及机身坡度综合位移施加图，此时在盾为盾构完全进入土层后，超挖及机身坡度综合位移施加图，此时在盾构机尾内部已完成拼装一环衬砌。构
44、机尾内部已完成拼装一环衬砌。整环衬砌安装完毕图图(b)为盾构开挖至第二环距离时，盾构机身相应位置施加位移为盾构开挖至第二环距离时，盾构机身相应位置施加位移 dr、2dr2.3.10 注浆材料时空变化特性模拟注浆材料时空变化特性模拟图图(d)为第一环衬砌安装完毕，进行同步注浆。此时，将注浆单元层位为第一环衬砌安装完毕，进行同步注浆。此时，将注浆单元层位置的机身单元置的机身单元“杀死杀死”，模拟盾尾前移，并解除该处机身与土体的接触。，模拟盾尾前移，并解除该处机身与土体的接触。考虑到注浆力学性质随时间和空间逐步变化，因此将注浆考虑成非线性考虑到注浆力学性质随时间和空间逐步变化，因此将注浆考虑成非线性
45、弹性材料，其弹性模量随时间变化。弹性材料，其弹性模量随时间变化。实际上，注浆材料在随时间变化过程中，要经历液态到固态的转变，通实际上，注浆材料在随时间变化过程中，要经历液态到固态的转变，通过在数值模型中改变注浆层单元弹性模量的方法，一定程度上反映了在液态过在数值模型中改变注浆层单元弹性模量的方法，一定程度上反映了在液态阶段注浆与土体的相互融合作用，产生较大变形；在固态阶段随着水合作用阶段注浆与土体的相互融合作用，产生较大变形；在固态阶段随着水合作用不断加强，注浆变形逐渐趋于稳定的过程。不断加强，注浆变形逐渐趋于稳定的过程。2.3.9 刀盘切削作用模拟刀盘切削作用模拟盾构刀盘扭矩的绝大部分贡献在
46、了刀盘对开挖面土体切削作用上，盾构刀盘扭矩的绝大部分贡献在了刀盘对开挖面土体切削作用上，反映在有限元模型中，可在开挖面上施加均布的周向剪应力反映在有限元模型中，可在开挖面上施加均布的周向剪应力 ，且剪，且剪 应应力与刀盘扭矩相关，以此来模拟刀盘切削土体运动对地表沉降及开挖面力与刀盘扭矩相关，以此来模拟刀盘切削土体运动对地表沉降及开挖面稳定性的影响。该剪应力对开挖面中心的合扭矩即为刀盘扭矩，其大小稳定性的影响。该剪应力对开挖面中心的合扭矩即为刀盘扭矩，其大小计算如下式计算如下式式中T 为刀盘扭矩；R 为开挖面直径如图如图 2-10 所示为有限元模型的端所示为有限元模型的端面图。在未开挖前，衬砌及
47、注浆层预设面图。在未开挖前，衬砌及注浆层预设单元属性为相应位置土层特性。随着开单元属性为相应位置土层特性。随着开挖的进行，逐渐挖的进行，逐渐“激活激活”预设衬砌与注浆预设衬砌与注浆单元层的衬砌特性与注浆特性。单元层的衬砌特性与注浆特性。2.3.13 有限元模型及边界荷载施加有限元模型及边界荷载施加如图如图 2-11 所示为有限元模型注所示为有限元模型注浆层的材料属性随着盾构掘进时间与浆层的材料属性随着盾构掘进时间与空间位置的不同，其力学性质也不同空间位置的不同，其力学性质也不同（图中以不同颜色表示）（图中以不同颜色表示）随着开挖面前移，随着开挖面前移，开挖面泥水压力开挖面泥水压力作用也随之前移
48、。在有限元模型中，以作用也随之前移。在有限元模型中，以施加施加单元表面正压力单元表面正压力模拟，如图模拟，如图 2-12 所所示。示。开挖面处刀盘对土体开挖面处刀盘对土体切削作用切削作用采用采用施加施加表面剪应力表面剪应力的方式模拟，如图的方式模拟，如图 2-13 所示。开挖面泥水压力与刀盘扭矩剪应所示。开挖面泥水压力与刀盘扭矩剪应力在施加到开挖面土体的同时，也同时力在施加到开挖面土体的同时，也同时反作用到盾构前部泥水舱面板上反作用到盾构前部泥水舱面板上在有限元模型中，每组推进在有限元模型中，每组推进千斤顶推力千斤顶推力以以面压力的方式施加在面压力的方式施加在衬砌单元衬砌单元上，同时又反作用到
49、相应上，同时又反作用到相应位置的机身单元上，如图位置的机身单元上，如图 2-14 在盾构系统中，在盾构系统中，拖车拖车紧随盾构机，在已铺成的衬砌上前行，对隧道形紧随盾构机，在已铺成的衬砌上前行，对隧道形成一种移动荷载。在有限元模型上，通过图成一种移动荷载。在有限元模型上，通过图 2-15 来反映这种空间分布的移来反映这种空间分布的移动荷载，即以动荷载，即以正压力形式作用正压力形式作用在衬砌单元表面，并随盾构前行而前行。在衬砌单元表面，并随盾构前行而前行。盾构盾构机身机身在有限元模型中以在有限元模型中以壳单元壳单元模拟，并在不同的盾构机身位置施模拟，并在不同的盾构机身位置施加一定位移来模拟机身坡
50、度，如图加一定位移来模拟机身坡度，如图 2-16 所示所示如图为某一掘进时刻，各施工因素在有限元模型中的反映。当隧道穿越透水性如图为某一掘进时刻，各施工因素在有限元模型中的反映。当隧道穿越透水性较强的地层时，应考虑地下水对隧道的作用力且该作用力应随着隧道埋深不同而变较强的地层时，应考虑地下水对隧道的作用力且该作用力应随着隧道埋深不同而变化。在有限元模型中，化。在有限元模型中，地下水压力地下水压力以以面压力荷载面压力荷载形式，在穿越透水性较强土层时施形式，在穿越透水性较强土层时施加到加到衬砌层单元表面衬砌层单元表面，并随埋深不同而变化，并随埋深不同而变化.三、隧道衬砌的模拟三、隧道衬砌的模拟实体