岩土工程数值计算PPT课件.ppt

1、*1岩土工程数值计算二零一一年三月二零一一年三月*2参考文献参考文献计算土力学计算土力学 上海科学技术出版社上海科学技术出版社 朱百里朱百里土工计算机分析土工计算机分析 中国建筑工业出版社中国建筑工业出版社 龚晓南龚晓南土工数值计算土工数值计算 中国铁道出版社中国铁道出版社 钱家欢钱家欢土工原理与计算土工原理与计算 水利电力出版社水利电力出版社 钱家欢钱家欢岩土工程有限元分析理论与应用岩土工程有限元分析理论与应用 科学科学出版社出版社 谢康和谢康和Finite element analysis in geotechnical engineering David M. Potts *3有限元法的

2、基本知识有限元法的基本知识 0.*41 1 自由度的编号自由度的编号 N = n + 1 READ(1,*) (IB3(I),I=1,NP3) READ(1,*) (LH(I),I=1,NP5) DO 60 I=1,NP3 IS=IB3(I) CALL HDIR(IS) DO 60 J=2,3 IP3(JP(1),J-1)=JP(J)60 CONTINUE DO 70 I=1,NP5 IS=LH(I)/10 II=LH(I)-IS*10 IP3(IS,3)=II70 CONTINUE N3=0 DO 80 I=1,NP DO 100 J=1,3 IF(IP3(I,J).EQ.0) GOTO 1

3、00 IP3(I,J)=IP3(I,J)+N3 N3=N3+1100 CONTINUE80 CONTINUE213 0 1节点号节点号 X 方向约束信息方向约束信息Y 方向约束信息方向约束信息121 0 节点号节点号 0 代表孔压代表孔压*52 2 半带宽半带宽 dmdm dm = b + 1 b b为相邻自由度编号的最大差值（对应节点在同一单元）为相邻自由度编号的最大差值（对应节点在同一单元） nnmmkkkkkkkkkkkkkkkkkK99887766555453524544353325221100dm*6 IJK=3 NN=12 DO 10 I=1,4 IVP=IVEN(I) DO 10

4、 J=1,IJK IJ=IJK*(I-1)+J ICN(IJ)=IP(IVP,J)10 CONTINUE DO 20 J=1,NN ICNP=ICN(J) DO 20 K=1,NN ICNQ=ICN(K) IF(ICNP.EQ.0.OR.ICNQ.EQ.0) GOTO 20 IF(IA3(ICNP).GT.(ICNP-ICNQ) GOTO 20 IA3(ICNP)=ICNP-ICNQ20 CONTINUEHALF BAND WIDTH -1 RELATED TO Ith FREE DEGREE*73 3 总刚的组装总刚的组装l 节点平衡法节点平衡法12345每一单元节点力与节点位移之间的关系为

5、每一单元节点力与节点位移之间的关系为 22222RFFFFe321333231232221131211321kkkkkkkkkFFF5214212213212211211kkkkkkF3212211213212211212kkkkkkF2213215213212211213kkkkkkF根据外荷载与节点力平衡条件根据外荷载与节点力平衡条件*8l 直接刚度法直接刚度法先把每个单元的单刚阶数扩大成总刚阶数，把单刚中按局部编码的先把每个单元的单刚阶数扩大成总刚阶数，把单刚中按局部编码的子块搬到总刚中相应的总体编码的位置中去，余下部分用零子块填充。子块搬到总刚中相应的总体编码的位置中去，余下部分用零子

6、块填充。123456 32100000000000000000000000000654321333231232221131211kkkkkkkkkK1 2 3 4 5 61 2 3 *94 4 总刚的存储总刚的存储l 一维变带宽存储一维变带宽存储由于整体刚度矩阵具有对称性、稀疏性和带状性，采用变带宽下三由于整体刚度矩阵具有对称性、稀疏性和带状性，采用变带宽下三角一维存储。角一维存储。887877686766565546454435343324232214131211kkkkkkkkkkkkkkkkkkkkk887877686766565546454435343324232214131211kk

7、kkkkkkkkkkkkkkkkkkk*10 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 CONTINUE将某一单元的将某一单元的KeKe集整到劲度矩阵集整到劲度矩阵KK中中*11l SKYLINE ( 轮廓线法轮

8、廓线法 )按列存储刚度矩阵按列存储刚度矩阵上三角区上三角区必要部分，按行存储必要部分，按行存储下三角区下三角区中必要部中必要部分。对刚度矩阵中出现少数分。对刚度矩阵中出现少数非常长的列非常长的列的情况下，存储要求不会剧烈增的情况下，存储要求不会剧烈增加，很容易利用向量点积例行程序。加，很容易利用向量点积例行程序。8878776867665856554846454438353433282423221814131211kkkkkkkkkkkkkkkkkkkkkkkkkkK11 k12 k22 k13 k23 k33 k14 k24 k34*12DO 355 K=1,4 NRCC=3*(K-1) N

9、R=NQ(NOD(K,IX)-1 DO 350 M=1,2+IFLOW 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+

10、NCN S(NN-ISHIFT)=S(NN-ISHIFT)+C1(NRCC,NCCC) 344 CONTINUE 345 CONTINUE ENDIF SL(NR)=SL(NR) + ZY(NRCC) 350 CONTINUE 355 CONTINUE(NCN, NR)(NR, NCN)*135 5 边界条件的引入边界条件的引入l 划划 0 置置 1 法法： 处理处理 ui = 0 约束约束l 乘大数法：乘大数法：处理处理 ui = R 约束约束将总刚相应的主对角元素改为将总刚相应的主对角元素改为1 1，将对应的行、列其它元素改为，将对应的行、列其它元素改为0 0。将荷载向量中相应的元素改为将

11、荷载向量中相应的元素改为0 0。 或预先将每个节点的方程编号，已知或预先将每个节点的方程编号，已知位移的节点不编号。位移的节点不编号。将总刚相应的主对角元素置一大数，将荷载向量中相应的元素改为将总刚相应的主对角元素置一大数，将荷载向量中相应的元素改为该大数乘该大数乘 R R。*14l 初等变换法：初等变换法：处理处理 两个变量有确定关系两个变量有确定关系如轴对称问题，可以取一夹角为扇形区计算。在斜对称如轴对称问题，可以取一夹角为扇形区计算。在斜对称ABAB边上的点边上的点只能沿只能沿ABAB上移动，所以点上移动，所以点C C上的点的位移有上的点的位移有使用初等变换法使用初等变换法消除一个相关方

12、程消除一个相关方程，若消除，若消除对应的方程，可以将对应的方程，可以将所在的行乘所在的行乘tgtg 加到加到 u u 所在的行，并将所在的行，并将所在的列也乘所在的列也乘tgtg加到加到 u u 所在列上去。所在列上去。tgucc*15JP(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)=K3*166 6 有限元法解题步骤有限元法解题

13、步骤1 建立计算网格建立计算网格2 设定计算相关参数（计算精度、控制等）设定计算相关参数（计算精度、控制等）3 边界条件和内约束，力、位移、孔压等边界条件和内约束，力、位移、孔压等4 选择单元类型、本构模型，并输入本构参数选择单元类型、本构模型，并输入本构参数5 计算结果提取、分析等计算结果提取、分析等*17岩土工程问题分析方法岩土工程问题分析方法 1.*181.1 1.1 岩土工程问题控制方程的建立岩土工程问题控制方程的建立0uzyxzxyxx 0vzyxzyyxy 0gwzyxzyzxz 1 1 土体平衡方程土体平衡方程2 2 土体本构方程土体本构方程 Ddzzppdxxppgzxwwzp

14、p*193 3 土体几何方程土体几何方程xuxyvyzwzxvyuxyywzvyzzuxwzx4 4 土体有效应力原理土体有效应力原理pxxpyypzz*205 5 孔隙流体平衡方程孔隙流体平衡方程0ukvgxpwxxw 0vkvgypwyyw 0gwkvgzpwwzzw 6 6 渗流连续方程渗流连续方程zwyvxutzvyvxvzyxdzzzzdzxxxgzxw*217 7 总控制方程总控制方程uxpzxwddyxvddzudyudxud 2551324412225522442211)()(vypzywddyxuddzvdyvdxvd 2551324412225522112244)()(wg

15、zpzyvddyzuddzwdywdxwd 2551325513223322552255)()(0wkzvkyukxzwyvxutggzpkzypkyxpkxzyxwwwzyx *221.2 1.2 岩土工程基本分析方法岩土工程基本分析方法1 1 总应力分析法及其控制方程总应力分析法及其控制方程总应力分析法与一般固体力学相同。从应用上讲，一般用于不考虑渗总应力分析法与一般固体力学相同。从应用上讲，一般用于不考虑渗流固结的情况，如饱和粗粒土地基、透水性土料组成的土坝路堤的应力和流固结的情况，如饱和粗粒土地基、透水性土料组成的土坝路堤的应力和变形分析以及饱和软粘土地基短期变形和稳定性分析。变形分析

16、以及饱和软粘土地基短期变形和稳定性分析。0)()(2551324412225522442211zxwddyxvddzudyudxud0)()(2551324412225522112244zywddyxuddzvdyvdxvd0)()(2551325513223322552255zyvddyzuddzwdywdxwd*232 2 有效应力分析法及其控制方程有效应力分析法及其控制方程在有效应力分析法中，土体的有效应力和孔压被严格区分，并将土骨在有效应力分析法中，土体的有效应力和孔压被严格区分，并将土骨架变形与孔隙水的渗透同步考虑。因此，有效应力分析法较能更真实地反架变形与孔隙水的渗透同步考虑。因此

17、，有效应力分析法较能更真实地反映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广。映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广。有效应力分析法尚需要有效应力原理和连续性方程。总应力法中只有有效应力分析法尚需要有效应力原理和连续性方程。总应力法中只有位移变量且仅与空间有关，而有效应力法中还有孔压变量，而且与空间和位移变量且仅与空间有关，而有效应力法中还有孔压变量，而且与空间和时间均有关。时间均有关。0)()(2551324412225522442211xpzxwddyxvddzudyudxud0)()(2551324412225522112244ypzywddyxu

18、ddzvdyvdxvd0)()(2551325513223322552255gzpzyvddyzuddzwdywdxwdzwyvxutggzpkypkxpkwwzyx222222*243 3 总应力分析法总应力分析法和和有效应力分析法有效应力分析法关系关系总应力法总应力法是是有效应力法有效应力法中当孔压中当孔压p=0p=0时的特殊形式。在有效应力分析时的特殊形式。在有效应力分析中如果采用与总应力分析相同的土工参数，并令孔压中如果采用与总应力分析相同的土工参数，并令孔压p=0p=0，所得结果即为，所得结果即为总应力分析结果。总应力分析结果。总应力分析总应力分析一般采用土体得一般采用土体得不排水不

19、排水指标，由此进行的是加荷瞬时或短指标，由此进行的是加荷瞬时或短期应力和变形分析。但也可采用土体的排水指标，此时进行的是期应力和变形分析。但也可采用土体的排水指标，此时进行的是最终或长最终或长期期应力和变形分析。当进行线弹性分析并采用排水指标，总应力分析得到应力和变形分析。当进行线弹性分析并采用排水指标，总应力分析得到的结果为有效应力分析的最终（孔压消散完毕、主固结完成）结果。的结果为有效应力分析的最终（孔压消散完毕、主固结完成）结果。*254 4 不排水孔压计算不排水孔压计算 f 000ffffpppDffD fDDD333300IIKDef有效应力原理有效应力原理其中其中*26fsffpK

20、npKn1vv)(ezyxefKKpsfeKnKnK11孔压增量孔压增量不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须通过有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须通过颗粒接触，但接触面积很小，导致体积变化也很小，如果忽略这种体积变颗粒接触，但接触面积很小，导致体积变化也很小，如果忽略这种体积变化。则总的体积变化为化。则总的体积变化为*27对于对于饱和土体饱和土体，Ks、Kf 都比土骨架模量大很多，在不考虑具都比土骨架模量大很多，在不考虑具体值时（此时准确值已

21、不重要），可假设体值时（此时准确值已不重要），可假设Kf=Ks,得得feKK nKKfe由于由于Ks比土骨架模量大很多，如孔隙流体压缩性较大，以致比土骨架模量大很多，如孔隙流体压缩性较大，以致 Ks kf，则，则对于对于排水分析排水分析，取，取 Ke=0，加载过程中孔压不变。，加载过程中孔压不变。*28)21 (AAu3)1 ()21 ()1 (A对于各向同性线弹性土体对于各向同性线弹性土体进行进行不排水分析时，必须设置不排水分析时，必须设置Ke。据经验，对于饱和土。据经验，对于饱和土体，只要体，只要ke足够大，土体对足够大，土体对ke的实际大小并不敏感。但的实际大小并不敏感。但Ke太太大时，

22、可能导致数值不稳定，即不排水泊松比接近大时，可能导致数值不稳定，即不排水泊松比接近0.5。专家建议专家建议设置设置 ke = Kskel， 在在100和和1000之间，之间，Kskel是土骨架的体积模量。是土骨架的体积模量。*29u=0.1=0.3100.45200.47931000.49460.497710000.49940.4998)21 (AAu3)1 ()21 ()1 (A*301.3 1.3 岩土工程问题的边界条件岩土工程问题的边界条件1 1 固结分析中的边界条件固结分析中的边界条件对节点位移和孔压已知的情况，可以有两种处理方法对节点位移和孔压已知的情况，可以有两种处理方法（1 1）

23、仍给以自由度编号，在解方程时把解得的位移值用已知值置换。仍给以自由度编号，在解方程时把解得的位移值用已知值置换。（2 2） 已知值对应的自由度均编为已知值对应的自由度均编为0 0，把代数方程区分为：，把代数方程区分为：212122211211bbxxAAAA式中，为已知的边界值。计算中只对式中，为已知的边界值。计算中只对x x1 1建立方程，即建立方程，即 2121111xAbxA*31Displacement + pore pressure DOFDisplacement DOF其中，其中，A12x2项在对各单元进行计算时即可求出。在多数情况下，已知项在对各单元进行计算时即可求出。在多数情况

24、下，已知的边界值为的边界值为0 0，上式简化为，上式简化为1111bxA*322 2 渗流问题的边界条件渗流问题的边界条件岩土工程种进行渗流计算的主要目的是求得渗流量和孔隙水头岩土工程种进行渗流计算的主要目的是求得渗流量和孔隙水头或水力坡降。渗流边界可分为两类，即已知水头边界和已知流量边界或水力坡降。渗流边界可分为两类，即已知水头边界和已知流量边界，分别可以表示如下：，分别可以表示如下：hh1qhnhk2在渗流计算中，对于在渗流计算中，对于自由面渗流问题自由面渗流问题，其，其边界条件需特殊处理边界条件需特殊处理。*33(1) (1) 有压渗流有压渗流有压渗流的计算域边界和边界条件是固定的，可以

25、很方便地应用有限有压渗流的计算域边界和边界条件是固定的，可以很方便地应用有限元法计算。如图，元法计算。如图，abab和和efef为第一类边界，为第一类边界，cdcd和和fafa为为第二类边界第二类边界，bcbc和和eded如如果离所讨论问题的渗流区域足够远时，既可当作果离所讨论问题的渗流区域足够远时，既可当作第一类边界第一类边界，也可当作第，也可当作第二类边界。设二类边界。设n ni i为计算域内部节点数，为计算域内部节点数，nb1nb1和和nb2nb2分别为第一类和第二类的分别为第一类和第二类的节点数，节点数，nt=ni+nb2，则渗流的有限元公式可表示为：，则渗流的有限元公式可表示为：31

26、33injjijfhkbcdeh0h=h0h=h1h1a0q0q0qh=h1h=h0不透水层或或f*34或考虑或考虑 n=nt+nb1，并去掉下标后改写为：，并去掉下标后改写为：11133nbjjijinjjijhkfhktf fi i为为i i节点的节点的净流出水量净流出水量，对已知入渗的边界点，对已知入渗的边界点fifi是给定的，对内部节点是给定的，对内部节点和不透水边界点均为零。和不透水边界点均为零。bcdeh0h=h0h=h1h1a0q0q0qh=h1h=h0不透水层或或f*35水头求出后，通过单元中某一截面的渗流量为：水头求出后，通过单元中某一截面的渗流量为：ppniiixzniii

27、zxhzNlkhxNlkQ11其中其中np为单元内节点数，为单元内节点数，lx 和和 lz为计算截面在为计算截面在x x轴和轴和z z轴的投影。轴的投影。bcdeh0h=h0h=h1h1a0q0q0qh=h1h=h0不透水层或或f*36(2)(2)自由面渗流自由面渗流在堤坝、渠道和地下水渗流计算中，自由面或浸润面事先是未知的，在堤坝、渠道和地下水渗流计算中，自由面或浸润面事先是未知的，因此计算域本身是未知的。因此计算域本身是未知的。早期的研究工作是用早期的研究工作是用试算法试算法不断改变自由面位置不断改变自由面位置，直到计算结果满足，直到计算结果满足如下自由面上的边界条件为止。如下自由面上的边

28、界条件为止。zh (1)(2)0nh*37然后用新坐标进行下一轮计算，直到自由面上所有节点满足然后用新坐标进行下一轮计算，直到自由面上所有节点满足 h-z h-z pfi 流量边界Pfbpfi 孔压边界*161盾构隧道施工过程分析盾构隧道施工过程分析 9.*162E D C B APost shield lossshield lossface lossDirection of driveA cutting head and hoodB shieldC tailD unground liningE grouted lining一、隧道施工的地面响应一、隧道施工的地面响应地层损失：地层损失： 由于

29、地层向开挖面移动，额外被开挖的超过隧道体积的土体，通常用隧道体积的百分比表示。*1631. Gap methodgapDmDl二、隧道开挖的模拟二、隧道开挖的模拟Rowe etal 1983Initial tunnelpositionFinal tunnelpositionThe gap=Dm-D1 at the crown,and is always fixed to zero at the invert.在有限元网格中预先预先定义空隙，代表预期的地层损失。通过改变空隙的尺寸，可以反映不同的隧道施工方式。空隙在隧道周围设置，因而确定了开挖前的土体边界。将隧道仰拱设在下卧土体上，在隧道拱顶规定

30、空隙参数。*164gapDmDlInitial tunnelpositionFinal tunnelpositionThe gap=Dm-D1 at the crown,and is always fixed to zero at the invert.分析时，去掉开挖面周边的剪力，并监测节点位移。当节点位移达到空分析时，去掉开挖面周边的剪力，并监测节点位移。当节点位移达到空隙合拢、土体与预先定义的衬砌位置接触时，在这些节点激活土体隙合拢、土体与预先定义的衬砌位置接触时，在这些节点激活土体衬砌相衬砌相互作用。互作用。*1652. Convergence-confinement methodin

31、itialLining installation/unloadingExcavation/unloading 0i 0dd 0*)1 (ddPanet and Guenot 1982设定衬砌施工前的卸载比例，因此体积损失是一个预测值。在隧道边界节点施加力矢量（1-）F0（ F0 相当于土体初始应力0 ）初始等于0，以后逐渐增加到1 以模拟开挖过程。 在规定值d时，安装衬砌，这时边界应力减小d 0 ，剩余应力部分新增了衬砌应力。*166Ee3. Progressive softening methodinitialExcavation, Lining installationand unload

32、ingStiffness reductionand unloading0EEe0EEeeEE0E0 0 *Swoboda 1979通过给土体刚度乘一降低因子降低因子，使掘进头里的土体软化。当开挖力作用于未来隧道边界时，土体出现明显软化。当隧道采用台阶法施工时，每级施工都可用上述方法。*1674. Volume loss methodVsVs=volume of settlement troughzVL = Ve - VtVe = excavated tunnel volumeVt = final tunnel volume类似于收敛限制法，但不规定衬砌施工前的卸载比例，而是规定开挖完成后的体积

33、损失。此法用于体积损失可以预期的给定的掘进法*168Forces imposedBy soil to be excavatedForces imposed ontunnel boundary due to excavationF-F 计算由被开挖土体施加于隧道边界的等效节点力F0，它被线性分配到开挖发生过程的n个增量步，即F=F0/n。 在开挖的每一增量步在开挖边界施加等值反向的力矢量-F。*169 监测每一步的体积损失，在一预想的体积损失增量步施工衬砌。 衬砌施工后，对于剩下的增量步，仍然在开挖周边施加载荷-F，因而在衬砌内引起初始应力。后续过程还有可能发生进一步的地层损失，这取决于衬砌的刚

34、度。因此，在小于预期地层损失的增量步安装衬砌，完成开挖后可以得到预期的地层损失。Forces imposedBy soil to be excavatedForces imposed ontunnel boundary due to excavationF-F*170开挖面释放荷载的计算开挖面释放荷载的计算在盾构隧道施工过程中，由于土压舱内土体压力不能完全平衡开在盾构隧道施工过程中，由于土压舱内土体压力不能完全平衡开挖面土压力，使得土体开挖面前方土体发生趋近或者远离盾构开挖面挖面土压力，使得土体开挖面前方土体发生趋近或者远离盾构开挖面的三维位移。的三维位移。现有的模拟方法是每向前开挖一步，通过

35、开挖面处的土体原位土现有的模拟方法是每向前开挖一步，通过开挖面处的土体原位土压力与土压舱支护压力的差值压力与土压舱支护压力的差值P 计算作用在开挖面结点上的等效节计算作用在开挖面结点上的等效节点力。这种方法的开挖面前方土体主要发生点力。这种方法的开挖面前方土体主要发生水平方向的位移水平方向的位移,而实际上，而实际上，由于开挖面卸荷引起的前方由于开挖面卸荷引起的前方土体移动是三维土体移动是三维的的 ，因此，该处理方法并，因此，该处理方法并不合理。不合理。*171本文采用基于本文采用基于 Mana 法法 的开挖释放荷载计算方法如图所示。将开的开挖释放荷载计算方法如图所示。将开挖面正前方一层土体单元

36、作为开挖卸荷单元，用该层开挖卸荷单元的收挖面正前方一层土体单元作为开挖卸荷单元，用该层开挖卸荷单元的收缩变形反映开挖引起的土体位移。缩变形反映开挖引起的土体位移。Mana 法的表达式为法的表达式为*172由于盾构施工中开挖面上应力由于盾构施工中开挖面上应力没有完全释放没有完全释放，因而不能像模拟基坑，因而不能像模拟基坑开挖那样，用上一步得到的单元应力直接计算新开挖面处的开挖释放荷开挖那样，用上一步得到的单元应力直接计算新开挖面处的开挖释放荷载，而应该用单元载，而应该用单元应力改变量应力改变量计算开挖释放荷载。计算开挖释放荷载。本文认为由于刀盘的搅动开挖面前方土体受到很大的扰动，可近似本文认为由

37、于刀盘的搅动开挖面前方土体受到很大的扰动，可近似认为土体由原来的的应力状态变为各向等压状态，并且压力等于土压舱认为土体由原来的的应力状态变为各向等压状态，并且压力等于土压舱支护压力支护压力 ps ，土体单元应力改变量，土体单元应力改变量同时将开挖卸荷单元的模量降低。同时将开挖卸荷单元的模量降低。x， y， z = ps- x，ps-y， ps-z*173分别运用直接由分别运用直接由P 计算等效节点力的方法和本文推荐方法，计算计算等效节点力的方法和本文推荐方法，计算得到纵断面上由开挖面卸荷引起的盾构周围土体位移场如图得到纵断面上由开挖面卸荷引起的盾构周围土体位移场如图 2 所示（计所示（计算中采

38、用的支护压力均小于开挖前土体初始应力），可以看出本文推荐算中采用的支护压力均小于开挖前土体初始应力），可以看出本文推荐的方法能够较好的模拟开挖面前方土体的三维移动。的方法能够较好的模拟开挖面前方土体的三维移动。*174刀盘超挖间隙和盾尾空隙的模拟刀盘超挖间隙和盾尾空隙的模拟为了使盾构能够顺利前行而减小盾构壳上的摩擦阻力，通常盾构的为了使盾构能够顺利前行而减小盾构壳上的摩擦阻力，通常盾构的刀盘外径要大于盾构壳的外径，从而在盾构壳外围产生一定厚度的间隙刀盘外径要大于盾构壳的外径，从而在盾构壳外围产生一定厚度的间隙 ，盾构机由于自重发生下沉到底部，在横断面上形成超挖间隙，随后，盾构机由于自重发生下沉

39、到底部，在横断面上形成超挖间隙，随后，周 围 土 体 由 于 填 充 超 挖 间 隙 而 产 生 指 向 盾 构 内 部 的 径 向 位 移周 围 土 体 由 于 填 充 超 挖 间 隙 而 产 生 指 向 盾 构 内 部 的 径 向 位 移*175当盾构前行时，由于不能做到完全同步注浆，而在盾尾产生建筑空当盾构前行时，由于不能做到完全同步注浆，而在盾尾产生建筑空隙，如图隙，如图 4 所示。图中所示。图中 t 为盾构壳的厚度，空隙所对应的衬砌管片由于为盾构壳的厚度，空隙所对应的衬砌管片由于受到前后衬砌螺栓的固定作用不能下沉，因而形成图示的盾尾建筑空隙。受到前后衬砌螺栓的固定作用不能下沉，因而形

40、成图示的盾尾建筑空隙。*176模拟上述刀盘超挖间隙和盾尾建筑空隙时，首先将相应位置的盾模拟上述刀盘超挖间隙和盾尾建筑空隙时，首先将相应位置的盾构壳单元和衬砌单元用刚度极小的软材料代替（活化减退），给外围构壳单元和衬砌单元用刚度极小的软材料代替（活化减退），给外围土体施加如图所示的已知径向位移，当土体接触到盾构壳或者衬砌外土体施加如图所示的已知径向位移，当土体接触到盾构壳或者衬砌外壁时，再将盾构壳和衬砌单元重新激活，使它们与土体共同作用壁时，再将盾构壳和衬砌单元重新激活，使它们与土体共同作用*177盾尾注浆的模拟采用对盾构壳外围的土体单元施加远离隧道中心盾尾注浆的模拟采用对盾构壳外围的土体单元施

41、加远离隧道中心的结点荷载，其数值由注浆压力计算得到的结点荷载，其数值由注浆压力计算得到盾尾注浆的模拟盾尾注浆的模拟F盾构法隧道施工的精细模拟，张海波，岩土力学*178刚度迁移法是将盾构前行看成刚度和载荷的迁移过程。图刚度迁移法是将盾构前行看成刚度和载荷的迁移过程。图2是盾是盾构前行时刚度迁移法示意图。在盾首和盾尾均设有预设单元，开挖面构前行时刚度迁移法示意图。在盾首和盾尾均设有预设单元，开挖面推进时，盾首逐渐深入，预设单元刚度逐渐增加推进时，盾首逐渐深入，预设单元刚度逐渐增加;盾尾逐渐脱出，预盾尾逐渐脱出，预设单元的刚度逐渐减少。当推进了一段行程后，盾首预设单元变为外设单元的刚度逐渐减少。当推

42、进了一段行程后，盾首预设单元变为外壳刚度，盾尾的变为空气单元，盾构前行了一个行程。壳刚度，盾尾的变为空气单元，盾构前行了一个行程。刚度迁移法原理刚度迁移法原理*179盾构前行的同时，盾构附属的其他结构包括载荷也必须前行。盾构前行的同时，盾构附属的其他结构包括载荷也必须前行。加强环迁移与预设单元前行相同，载荷的迁移按加强环的刚度迁移加强环迁移与预设单元前行相同，载荷的迁移按加强环的刚度迁移来分配。当前加强环由前位置迁移到后位置时，载荷也随着迁移。来分配。当前加强环由前位置迁移到后位置时，载荷也随着迁移。最后，前盾构前行时，加强环、其上安装的设备及其载荷加强环成最后，前盾构前行时，加强环、其上安装

43、的设备及其载荷加强环成为空气单元时，其上的载荷也为零。为空气单元时，其上的载荷也为零。*180混合刚度接触单元混合刚度接触单元盾构推进过程的有限元模拟必须解决以下问题盾构推进过程的有限元模拟必须解决以下问题:(l) 盾构运动。盾构在推进力作用下向前运动，但是单元网格是业已划盾构运动。盾构在推进力作用下向前运动，但是单元网格是业已划好的，仅用盾构与土体之间的接触单元错动来模拟盾构前行，与盾构好的，仅用盾构与土体之间的接触单元错动来模拟盾构前行，与盾构前行不符，见图前行不符，见图1。刚度迁移法*181(2) 运动接触点对。盾构前行产生很大的相对位移，盾构对运动接触点对。盾构前行产生很大的相对位移，

44、盾构对B点产生支承，点产生支承，支承面支承面BB，AB成为自由面，失去盾构的支承，成为自由面，失去盾构的支承，A不再与不再与A点形成结点对。点形成结点对。(3)力学性质改变。盾构与土体间的滑动过程是有摩擦滑动，滑动在力学性质改变。盾构与土体间的滑动过程是有摩擦滑动，滑动在AB间产生，但当间产生，但当A了已经脱出盾尾时，此间不会产生摩擦。了已经脱出盾尾时，此间不会产生摩擦。*182随着盾构迁移，盾首和盾尾与土体的接触点在移动。以盾尾为例，随着盾构迁移，盾首和盾尾与土体的接触点在移动。以盾尾为例，接触点对接触点对AB运动的结果是最后一层，盾尾单元一部分是外壳刚度，运动的结果是最后一层，盾尾单元一部

45、分是外壳刚度，另一部分则是空气单元刚度，两种不同的材料填充在一个单元内，见另一部分则是空气单元刚度，两种不同的材料填充在一个单元内，见图。这样的特殊单元是不能直接用图。这样的特殊单元是不能直接用Goodman单元来模拟的，必须加单元来模拟的，必须加以修正。以修正。盾尾两种介质的混合单元*183在盾构迁移过程中，盾首和盾尾单元及其接触单元内一部分由土在盾构迁移过程中，盾首和盾尾单元及其接触单元内一部分由土体材料变为外壳材料，或由外壳材料变为空气材料，见图体材料变为外壳材料，或由外壳材料变为空气材料，见图3。单元内两。单元内两种材料影响到应力一应变矩阵，即单元刚度。假设己知单元结点的刚种材料影响到

46、应力一应变矩阵，即单元刚度。假设己知单元结点的刚度，那么单元内任意一点刚度可以用加权形函数插值求得。度，那么单元内任意一点刚度可以用加权形函数插值求得。(l) 空间空间Goodman单元内两种材料下任一点权刚度单元内两种材料下任一点权刚度式中：式中： i=Ai/A为为Ei面积权系数，反映盾尾位置对单元的影响面积权系数，反映盾尾位置对单元的影响*184(2) 盾构外壳单元内两种材料的任一点权弹性常数盾构外壳单元内两种材料的任一点权弹性常数式中：式中： i=Ai/A为为Ei面积权系数，反映盾尾位置对单元的影响面积权系数，反映盾尾位置对单元的影响*185从盾构法隧道施工过程的整体来看，实际上是一种从

47、盾构法隧道施工过程的整体来看，实际上是一种稳态掘进稳态掘进的过程，的过程，在每一时刻，盾构系统与周围土层均处于在每一时刻，盾构系统与周围土层均处于平衡状态平衡状态。然而，从。然而，从局部局部来看，来看，特别是盾构切削土体这一过程，是一种复杂地特别是盾构切削土体这一过程，是一种复杂地动态相互作用动态相互作用过程。因此，过程。因此，要研究盾构机的掘进性能，必然要考虑到盾构掘进过程中的宏观平衡特要研究盾构机的掘进性能，必然要考虑到盾构掘进过程中的宏观平衡特性以及局部动态特性。性以及局部动态特性。目前，国内外对盾构目前，国内外对盾构刀盘切削土体刀盘切削土体过程的数值模拟研究非常少，过程的数值模拟研究非

48、常少，这主要由于：一方面，要实现刀盘对土体的切削模拟，必须能够解决这主要由于：一方面，要实现刀盘对土体的切削模拟，必须能够解决刀盘与土体的刀盘与土体的接触行为接触行为模拟，处理由土体大变形、破坏等带来的数值模拟，处理由土体大变形、破坏等带来的数值稳定性；另一方面，刀盘对土体切削过程，实际上是一个稳定性；另一方面，刀盘对土体切削过程，实际上是一个 3 尺度问题尺度问题，即，刀齿切削土体过程为厘米级别，刀盘与土体相互作用过程为米级即，刀齿切削土体过程为厘米级别，刀盘与土体相互作用过程为米级别，整个地层环境为百米级别以上。因此，要对这一个多尺度问题进别，整个地层环境为百米级别以上。因此，要对这一个多

49、尺度问题进行三维动态模拟，既需要合适的算法支持，更需要强大的计算能力。行三维动态模拟，既需要合适的算法支持，更需要强大的计算能力。*186盾构掘进离散方法沈建奇，盾构掘进过程数值模拟方法研，上海交通大学博士学位论文，2009泥水盾构掘进过程的逐步离散，是将盾构掘进的泥水盾构掘进过程的逐步离散，是将盾构掘进的连续连续稳态过程稳态过程离散离散成以一环衬砌距离为成以一环衬砌距离为掘进步长掘进步长的的逐步逐步稳态掘进过程稳态掘进过程。将泥水盾构每掘进。将泥水盾构每掘进一环衬砌距离，即认为盾构一环衬砌距离，即认为盾构-土层土层-构筑物系统处于构筑物系统处于外部和内部平衡状态外部和内部平衡状态，将这一平衡

50、状态作为一载荷步进行迭代求解。将这一平衡状态作为一载荷步进行迭代求解。在盾构掘进过程中，相应的土体单元被在盾构掘进过程中，相应的土体单元被“杀死杀死”，同时相应位置，同时相应位置的衬砌、注浆、机身单元被的衬砌、注浆、机身单元被“激活激活”，相应的荷载条件、边界条件被施，相应的荷载条件、边界条件被施加，并在每一掘进步内，计算盾构加，并在每一掘进步内，计算盾构-土体土体-构筑物系统的应力、应变，最构筑物系统的应力、应变，最后得到以掘进距离为横坐标的位移、变形、应力、姿态等重要数据变化后得到以掘进距离为横坐标的位移、变形、应力、姿态等重要数据变化曲线。曲线。*187下图表达了盾构机头进入土层，然后连