岩土工程数值计算

1、岩土工程数值计算,参考文献,计算土力学 上海科学技术出版社 朱百里,土工计算机分析 中国建筑工业出版社 龚晓南,土工数值计算 中国铁道出版社 钱家欢,土工原理与计算 水利电力出版社 钱家欢,岩土工程有限元分析理论与应用 科学出版社 谢康和,Finite element analysis in geotechnical engineering David M. Potts,有限元法的基本知识,0,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

2、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)=II 70 CONTINUE N3=0 DO 80 I=1,NP DO 100 J=1,3 IF(IP3(I,J).EQ.0) GOTO 100 IP3(I,J)=IP3(I,J)+N3 N3=N3+1 100 CONTINUE 80 CONTINUE,213 0 1,节点号 X 方向约束信息 Y 方向约束信息,121 0,节点号 0 代表孔压,2 半带宽 dm,dm = b + 1,b

3、为相邻自由度编号的最大差值（对应节点在同一单元,dm,IJK=3 NN=12 DO 10 I=1,4 IVP=IVEN(I) DO 10 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-ICNQ 20 CONTINUE,HALF BAND WIDTH -1 RELATE

4、D TO Ith FREE DEGREE,3 总刚的组装,节点平衡法,1,2,3,4,5,每一单元节点力与节点位移之间的关系为,根据外荷载与节点力平衡条件,直接刚度法,先把每个单元的单刚阶数扩大成总刚阶数，把单刚中按局部编码的子块搬到总刚中相应的总体编码的位置中去，余下部分用零子块填充,1,2,3,4,5,6,1 2 3 4 5 6,1 2 3,4 总刚的存储,一维变带宽存储,由于整体刚度矩阵具有对称性、稀疏性和带状性，采用变带宽下三角一维存储,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

5、 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,将某一单元的Ke集整到劲度矩阵K中,SKYLINE ( 轮廓线法,按列存储刚度矩阵上三角区必要部分，按行存储下三角区中必要部分。对刚度矩阵中出现少数非常长的列的情况下，存储要求不会剧烈增加，很容易利用向量点积例行程序,K11 k12 k22 k13 k23 k33 k14 k24 k34,DO 355

6、 K=1,4 NRCC=3*(K-1) NR=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)

7、NN=LOCC(NR)-LENNC-NR+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,5 边界条件的引入,划 0 置 1 法： 处理 ui = 0 约束,乘大数法：处理 ui = R 约束,将总刚相应的主对角元素改为1，将对应的行、列其它元素改为0。将荷载向量中相应的元素改为0。 或预先将每个节点的方程编号，已知位移的节点不编号,将总刚相应的主对角元素置一

8、大数，将荷载向量中相应的元素改为该大数乘 R,初等变换法：处理 两个变量有确定关系,如轴对称问题，可以取一夹角为扇形区计算。在斜对称AB边上的点只能沿AB上移动，所以点C上的点的位移有,使用初等变换法消除一个相关方程，若消除对应的方程，可以将所在的行乘tg 加到 u 所在的行，并将所在的列也乘tg加到 u 所在列上去,JP(1): NODE NUMBER JP(2): X-DIRECTION RESTRAIN INFORMATION JP(3): Y-DIRECTION RESTRAIN INFORMATION,JP(1)=I3/100 K3=I3-JP(1)*100 JP(2)=K3/10

9、K3=K3-JP(2)*10 JP(3)=K3,6 有限元法解题步骤,1 建立计算网格 2 设定计算相关参数（计算精度、控制等） 3 边界条件和内约束，力、位移、孔压等 4 选择单元类型、本构模型，并输入本构参数 5 计算结果提取、分析等,岩土工程问题分析方法,1,1.1 岩土工程问题控制方程的建立,1 土体平衡方程,2 土体本构方程,z,x,w,wz,p,p,3 土体几何方程,4 土体有效应力原理,5 孔隙流体平衡方程,6 渗流连续方程,z,x,w,7 总控制方程,1.2 岩土工程基本分析方法,1 总应力分析法及其控制方程,总应力分析法与一般固体力学相同。从应用上讲，一般用于不考虑渗流固结的

10、情况，如饱和粗粒土地基、透水性土料组成的土坝路堤的应力和变形分析以及饱和软粘土地基短期变形和稳定性分析,2 有效应力分析法及其控制方程,在有效应力分析法中，土体的有效应力和孔压被严格区分，并将土骨架变形与孔隙水的渗透同步考虑。因此，有效应力分析法较能更真实地反映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广,有效应力分析法尚需要有效应力原理和连续性方程。总应力法中只有位移变量且仅与空间有关，而有效应力法中还有孔压变量，而且与空间和时间均有关,3 总应力分析法和有效应力分析法关系,总应力法是有效应力法中当孔压p=0时的特殊形式。在有效应力分析中如果采用与总应力分析相同的土工参数，并

11、令孔压p=0，所得结果即为总应力分析结果,总应力分析一般采用土体得不排水指标，由此进行的是加荷瞬时或短期应力和变形分析。但也可采用土体的排水指标，此时进行的是最终或长期应力和变形分析。当进行线弹性分析并采用排水指标，总应力分析得到的结果为有效应力分析的最终（孔压消散完毕、主固结完成）结果,4 不排水孔压计算,有效应力原理,其中,孔压增量不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须通过颗粒接触，但接触面积很小，导致体积变化也很小，如果忽略这种体积变化。则总的体积变化为,对于饱和土体，Ks、Kf 都比土骨架模量大很多，在不考虑具体

12、值时（此时准确值已不重要），可假设Kf=Ks,得,由于Ks比土骨架模量大很多，如孔隙流体压缩性较大，以致 Ks kf，则,对于排水分析，取 Ke=0，加载过程中孔压不变,对于各向同性线弹性土体,进行不排水分析时，必须设置Ke。据经验，对于饱和土体，只要ke足够大，土体对ke的实际大小并不敏感。但Ke太大时，可能导致数值不稳定，即不排水泊松比接近0.5,专家建议设置 ke = Kskel， 在100和1000之间，Kskel是土骨架的体积模量,1.3 岩土工程问题的边界条件,1 固结分析中的边界条件,对节点位移和孔压已知的情况，可以有两种处理方法 （1）仍给以自由度编号，在解方程时把解得的位移值

13、用已知值置换。 （2） 已知值对应的自由度均编为0，把代数方程区分为,式中，为已知的边界值。计算中只对x1建立方程，即,Displacement + pore pressure DOF,Displacement DOF,其中，A12x2项在对各单元进行计算时即可求出。在多数情况下，已知的边界值为0，上式简化为,2 渗流问题的边界条件,岩土工程种进行渗流计算的主要目的是求得渗流量和孔隙水头或水力坡降。渗流边界可分为两类，即已知水头边界和已知流量边界，分别可以表示如下,在渗流计算中，对于自由面渗流问题，其边界条件需特殊处理,1) 有压渗流,有压渗流的计算域边界和边界条件是固定的，可以很方便地应用有

14、限元法计算。如图，ab和ef为第一类边界，cd和fa为第二类边界，bc和ed如果离所讨论问题的渗流区域足够远时，既可当作第一类边界，也可当作第二类边界。设ni为计算域内部节点数，nb1和nb2分别为第一类和第二类的节点数，nt=ni+nb2，则渗流的有限元公式可表示为,f,或考虑 n=nt+nb1，并去掉下标后改写为,fi为i节点的净流出水量，对已知入渗的边界点fi是给定的，对内部节点和不透水边界点均为零,f,水头求出后，通过单元中某一截面的渗流量为,其中np为单元内节点数，lx 和 lz为计算截面在x轴和z轴的投影,f,2)自由面渗流,在堤坝、渠道和地下水渗流计算中，自由面或浸润面事先是未知

15、的，因此计算域本身是未知的,早期的研究工作是用试算法不断改变自由面位置，直到计算结果满足如下自由面上的边界条件为止,1,2,然后用新坐标进行下一轮计算，直到自由面上所有节点满足 h-z 精度要求，可取最大水头差的1,计算中把自由面当作不透水面，所以式（2）自动满足，而式（1）一般不满足。 如果某一节点在第i次计算中假定的节点高程为zi，求得的节点水头为hi，不满足式（1），则利用下式计算新的坐标,具体计算时把渗流域分成非饱和的负压区、饱和的正压区和过渡区三部分，并把渗透系数随孔压的变化简化成三段折线，假定负压区的渗透系数很小，等于饱和区的1，过渡区内渗透系数随孔压p =g(h-z)直线变化。此

16、时下图边界条件为,h1,h0,g f,a,e,d,b c,ab 和 ag,bc 和 fg,cd,de 和 ef,另一种方法是把自由面以上的非饱和区和自由面以下的饱和区同时计算。计算所得的零孔隙压力线即为满足h=z条件的自由面,3) 不稳定渗流,随时间变化的渗流称为不稳定渗流。随动力荷载迅速变化的渗流问题需结合土骨架的变形同时计算。水库中水位降落时的渗透力将对土坝上游边坡造成不利影响，是不稳定渗流计算的主要对象,不考虑水体的压缩性，不稳定渗流控制方程仍如前所述，只是边界条件随时间变化。特别在自由面渗流中，自由面也是变化的，此时自由面除应满足 h=z 外，还应满足下列边界条件,式中ne为有效孔隙率

17、，相当于单位土体孔隙中自由水的含量,许多文献将上式表示为,此式的物理含义是，进入自由面的水量恰为自由面变动区的土体（图中阴影区）内所含的水量,ne 为给水度,把 写成差分形式： ，h0 为上一时段的水头值，并代入有限元插值公式 ，则某一节点i的渗出量为,nf为自由面节点数， 。如此，前述有限元公式可以改写为,如此，前述有限元公式 可以改写为,nf2 = nf - nf1，nf1 为自由面上水头已知节点数，nt中仍包括自由面节点nf2，nb1为包括nf1在内的所有水头已知节点。由于自由面上还需满足h=z的条件，所以对每一时段t还需进行自由面位置的迭代计算,f,真空堆载预压分析,2,2.1 施工工

18、艺,先打袋装砂井（直径约7cm，间距约1m），上面铺50cm的砂垫层，其中埋有真空滤管及与其连接的主管。然后在加固区四周开挖深约1m的沟槽，铺上塑料薄膜，将其周边埋入沟中，并将主管伸出薄膜与真空设备相连,密封沟,抽真空系统,主管,支 滤 管,1 2 3 4,1 2 3 4,2.2 控制方程,在平面上用有限单元法离散化，时间上用差分法分段以后，BIOT固结理论可以表示为,其中 、 和 为时段内j节点的水平、垂直位移和孔压增量， 、 为,i 结点的水平、垂直向荷载增量，pj0为上一时段末j结点的孔隙压力，为差分常数，此处用2/3。 等为方程式的系数，反映j结点对i结点的影响。其中 、 和 、 与选

19、用的单元形状有关， 除单元形状外，还与土的渗透系数有关，而 、 和 除单元形状外和，还与土骨架的变形性质有关,作为一个边值问题，土体内部发生固结变形过程是边界条件改变的结果。边界条件可以分为荷载边界、位移边界和水头边界、水量边界几种。真空预压与加荷预压只是在边界条件上有所不同，而控制方程式与解题方法是完全一样的,在加荷压缩问题中，只有边界 荷载（应力控制）或边界位移（应变控制）发生 变化，水流 边界条件 不变； 在真空抽水问题中，则只有水流边界条件发生变化，即真空作用点的水头下降，从而形成由内向外的水力坡降。显然，两者同时变化的问题（真空联合堆载）也可同样求解,说 明,一般固体力学问题的属性通

20、常仅指应变或应力属性，而固结问题的属性则包括应变和渗流属性。在固体力学中，平面问题指应变或应力是二维的问题，但对固结而言，当变形和渗流均为二维时才属于平面问题,2.3 固结问题属性,在一般的实际荷载作用下，天然软土地基或人工处理后地基固结问题通常属于变形和渗流均为三维的空间问题。但在长条形荷载作用下，天然软土地基固结问题属于平面变形、平面渗流固结问题（即平面问题），而经人工处理后的地基固结则往往属于平面变形、空间渗流固结问题,为了编程和计算方便，常将砂井固结问题作为一般的平面问题（平面变形、平面渗流）分析。将砂井简化为沿荷载长度方向连续（垂直于计算断面）的砂墙，并通过改变砂井和土体的渗透系数等

21、来使原砂井地基与简化后的砂墙地基等效,2.4 砂井（排水板）的处理,第种方法: 将竖向排水通道视为渗透性很强的介质,由于砂井的间距很小， 而为了得到合理的孔隙压力分布，两井之间至少应布置三排结点 这样势必使结点数成倍增加，增加计算工作量。因此还需要放大砂井的间距。 但这两种变换应保证变换前后主要基本量(如固结度或同一深度处的平均孔压)保持不变的前提下进行，这可以通过在平面应变有限元计算时调整砂井地基的实际渗透系数来实现,把土的水平向渗透系数按井距放大倍数的平方放大， 以保持水平向相对的渗径长度不变,这样的处理方法虽然比较粗糙 ，但经验表明 计算结果大体上还是能反映实际情况的。（,根据砂墙地基双

22、向应变、双向渗流等应变固结理论解与巴隆轴对称固结理论解的比较 ，在固结度或平均孔压不变的条件下，砂井地基与砂墙地基的等效可通过调整渗透系数得到 ，其等效计算公式为,Kxp = Dx kra,Kzp = Dz kza,Dx、Dz为水平向和竖向渗透系数的调整系数,kxp 、kzp分别为砂墙地基的水平和垂直渗透系数,kra 、kza 分别为砂井地基的水平和垂直渗透系数,实际计算时，用塑料排水板代替砂井，塑料排水板的当量直径采用下式 计算,式中,Dp 为塑料排水板当量直径,cm;为换算系数，无实验资料时可取0.751. 00; b为塑料排水板宽度; 为塑料排水板厚度,塑料排水板有效排水区直径按下式 计

23、算,式中, de为塑料排水板有效区直径;S 为塑料排水板间距,de = 1.05S,排水板与砂井的等效,在进行有限元计算时，根据实测膜下真空度值，将密封膜覆盖的地面边界单元结点的孔压设为负的稳定孔压（膜下真空度）。 在砂墙的位置，沿纵向与土接触的边界单元结点的孔压根据实测真空度调节为稳定负的孔压，从地面向下由大到小按线形变化设置，砂墙位置的底部负孔压为零（这与观测的孔隙水压力基本相符合) 。砂墙位置土的渗透系数及各种参数则跟周围土层相同,第 种方法: 将竖向排水通道视为负压边界条件,岑仰润.真空预压加固地基的试验及理论研究D.浙江大学 , 2003,设砂垫层和砂井中所有结点的孔隙水应力为常量，

24、砂垫层以外地表处的孔隙水应力为零，其他边界的孔隙水应力皆未知。设地表面自由变形,S1.透水边界 S2.负压边界 S3.自由边界 S4.不透水边界 S5.对称边界,由于竖向排水体井阻与涂抹作用的存在，抽真空作用在竖向排水体中形成的负压分布是复杂的，一些学者提出的真空预压加固地基时竖向排水体中的负压分布模式相差比较大,根据真空-堆载联合预压法的加固机理可知，真空荷载主要是作用在地表，并经由塑料排水板向下传递，使土体中的孔压与地表及排水板中的孔压存在压差而渗出，从而发生固结。 由于测量手段的局限，塑料排水板中的孔压变化规律尚不完全清楚，但是由于井阻等因素，塑料排水板中孔压的增长没有地表真空度那么迅速

25、和明显。由于塑料排水板中的孔压变化是由地表真空荷载变化引起，因此， 本文将地表的真空荷载作为已知边界条件，而将砂井中的孔压作为未知量处理。即将加固区地表的孔压变化值作为边界条件，堆载荷载作为一般荷载施加在加固区地表,第 III 种方法: 砂井孔压未知 (三维,孔压边界取为：砂垫层中所有结点的孔隙水压力为负的真空压力(80 kPa)，砂垫层以外的地基表面孔隙水压力为0，其它边界的孔压未知。 这是因为真空预压时，首先降低密封膜下砂垫层中的孔隙水压力，形成膜下真空度，并通过排水通道向下传递，使得土体内部各点与排水通道及砂垫层中形成压差，从而，发生由土中向边界渗流，而在固结过程中，竖向排水通道的孔压是

26、随时间变化的，孔隙水压力在加固过程中保持不变的只有砂垫层和砂垫层以外的地基表面,计算真空预压时，让加固区表面各点的孔压按线性从 0 减少到80 kPa，然后保持不变；在堆载时，将填土折算成等效结点荷载，按实际加载曲线施加上去。真空卸载时，让加固区表面各点的孔压按线性从80 kPa增加到0,竖向边界面： 水平位移约束、竖向位移无约束，不透水边界 基 底： 竖向、水平向位移约束，不透水边界。 地 面： 加固区和非加固区位移自由，非加固区孔压为0，加固区表面在抽真空时孔压大小为真空度，没有抽真空时该处的孔压也为0。真空度的变化和堆载的加载曲线，如图所示,具体而言，其约束条件为,假定在进行真空预压前，

27、各单元结点孔隙水应力皆为大气压力。为解方程的方便，我们以下都取相对压力，即设 pa = 0。于是有,式中0 非边界结点的孔隙水应力初始值。 设在抽气的瞬时，砂井及砂垫层中立即达到相同的负压值，且保持不变。 假定结点初始位移为零，即 0 = 0 式中0初始结点位移,0 = 0,2.5 初始条件,1 真空荷载处理,在砂垫层中抽真空时，负压作用在孔隙水上，使土中水部分流出土体，从而孔隙水应力降低，土骨架应力增加，而结点荷载始终为零，所以荷载向量 R=0 在解有限单元公式的方程组时，将负压边界(砂垫层、砂井）上的结点孔隙水应力的已知值乘以总系数矩阵中相应系数，移到方程右端作为荷载向量,2.6 荷载处理

28、,2 荷堆处理,对填土荷载的处理一般有2 种方法： 1）把填土当成新增单元参加计算； 2）将填土折算成等效结点荷载进行计算。 前者能够模拟填土与地基的相互作用，计算填土的应力和应变，分析粘性填土的变形和开裂问题具有明显的优势，但不能模拟线性加载过程；而后者可以方便地模拟实际的加载过程。 本文中采用第二种方式，填土荷载按施工设计分10 级施加，计算前将加载区域都设置为null（空单元），随后分别按时间顺序给所需加载的土层赋值，从而实现分级加载，计算模型也相应地分10 步完成,真空-堆载联合预压加固路堤FLAC数值分析,真空预压加固地基，均在地基中设置竖向排水体，竖向排水体的存在影响真空预压加固地

29、基的效果和地基的固结速度。竖向排水体一般采用袋装砂井和塑料排水板，早期多采用袋装砂井，如今基本上均采用塑料排水板。 当采用袋装砂井时，赋予砂井单元砂料的渗透系数，即己经考虑了井阻作用，若要再考虑涂抹作用，则将砂井单元周围土体划分为涂抹单元，涂抹区土体的强度和渗透系数适当折减以考虑涂抹作用。当采用塑料排水板时，可以将塑料排水板参照上述方法进行处理。从实际分析情况看，这样往往会过高估计负压沿竖向排水体的传递能力,2.7 竖向排水体的界定,考虑到竖向排水体井阻与涂抹作用的存在，抽真空作用在竖向排水体中形成一定的负压分布，也可认为地基是在竖向排水体中负压作用下固结的，因此也可以预先设定竖向排水体中负压

30、分布。从实际应用看，正确预先估计竖向排水体中负压分布模式是分析的关键,竖向排水体的界定,参考文献,2 软土地基真空排水预压的固结变形分析，岩土工程学报，1986，NO.3,1 用真空加固软土地基的机制与计算方法，岩土工程学报，1986，NO.2,3 刘汉龙，真空-堆载预压处理高速公路软基的有限元计算，岩土力学， 2003，No.6,4 岑仰润，真空预压加固地基的试验及理论研究D.浙江大学 , 2003,5 真空-堆载联合预压加固路堤FLAC数值分析,6 雷鸣，真空预压数值分析中竖向排水通道的讨论，铁道学报，2009,土与结构接触面单元,3,一、概 述,对于任何土和结构相互作用情况，结构与土之间

31、可能会发生相对运动,使用接触面单元，能模拟土与结构的接触边界条件。其优点是可以改变接触面的本构特性、允许土与结构之间的相对运动即滑动和分离。已有很多方法用于模拟土与结构接触面的非连续特性,当使用连续单元时，由于变形协调，会禁止土与结构接触面的相对运动如图，有限元法的节点协调条件将约束相邻的结构与土单元使其一起运动，这显然与实际不符,薄的连续单元， 标准的本构关系， Pande and Sharm 1979, Griffiths1985,2. 联动单元， 仅考虑相对节点之间的连接，通常采用独立的弹簧连接 Hermann 1978, Frank et al 1982,3. 接触面或节理单元， 零厚

32、度或有限厚度， Goodman et al 1968,，Ghaboussi et al 1973, Carol and Alonso 1983, Wilson 1977, Desai et al 1984, Beer 1985,4. 混合法,上述接触面的处理方法中，以零厚度接触面单元（Goodman 单元）应用最为普遍，后又被帝国理工学院Day1990予以发展,土和结构分开模拟，通过约束方程联系，使其在接触面上保持位移与力的协调性。 Francavilla and Zienkiewicz 1975, Katon 1983, Lai and booker 1989,二、无厚度接触面单元,Good

33、man等人提出一种接触面单元如图，由两片长度为L的接触面12和34组成。两片接触面间假想为无数微小的法向的和切向的弹簧连接。在受力前两接触面完全吻合，即单元没有厚度只有长度，是一维单元。一个接触面单元共4个节点如图,1 2,4 3,顶,底,z,x,1 2,4 3,顶,底,z,x,取线性位移模式，可将接触面上任一点的位移用节点位移表示,u3,v3,接触面单元内各点的相对位移为,式中,令,由虚位移原理可得,劲度系数（应力与相对位移的关系）取值,对于法向劲度系数，当接触面受压时，为了模拟两边二维单元不会在接触面处重叠，应取一很大的数值，如 kn = 108kN/m3，可使相互嵌入的相对位移小到可略去

34、。当接触面受拉时，因为接触面上不能承受拉应力，则令kn为很小的值，如取kn = 10kN/m3，以使算出的拉应力可忽略,拉 开,重 叠,对于切向劲度系数，可由直剪试验确定。Clough and Duncan摩擦试验表明，剪应力与相对剪切位移近似于双曲线关系，相应的表达式为,Clough and Duncan把初始剪切劲度和极限剪应力值与n的关系表示为,切线剪切劲度系数为,ksi,u,考虑孔压的单元刚度矩阵,K3,3=1 K6,6=1 K9,9=1 K12,12=1,K3,12 = -1 K6,9 = -1 K9,6 = -1 K12,3= -1,推求接触面单元的劲度矩阵Ke c SUBROUT

35、INE FJK(IE,JA) COMMON /PMEE/PME(19)/XZEE/XZE(4,2)/DDJJ/DJT(80,10) COMMON /EEKK/EK(12,12),ER(12) DO 40 I=1,12 DO 40 J=1,12 40 EK(I,J)=0.00 PA=98.10 G=1.00-DJT(IE,4)*PME(4) F=DJT(IE,10) IF(F.LT.0.50*PA) F=0.50*PA DJT(IE,1)=G*G*PME(5)*(F/PA)*PME(6)*0.10*PA IF(DJT(IE,4).GT.0.99990) DJT(IE,1)=0.10*PA DJT

36、(IE,2)=PME(7)*0.10*PA IF(DJT(IE,10).LE.0.1E-4*PA) DJT(IE,2)=0.10*PA SN=SIN(DJT(IE,3) CS=COS(DJT(IE,3) P=XZE(1,1)-XZE(2,1) Q=XZE(1,2)-XZE(2,2) A=SQRT(P*P+Q*Q) DO 60 M=1,4 DO 60 N=1,4 AKA=1.00,IF(M.GT.2.AND.N.LE.2.OR.M.LE.2.AND.N.GT.2) AKA=-1.00 AKB=1.00 IF(M.EQ.N.OR.(M+N).EQ.5) AKB=2.00 G=AKA*AKB*A/6

37、.00 I=JA*M+2-JA J=JA*N+2-JA EK(I-1,J)=G*(DJT(IE,1)-DJT(IE,2)*SN*CS EK(I,J-1)=EK(I-1,J) EK(I-1,J-1)=G*(DJT(IE,1)*CS*CS+DJT(IE,2)*SN*SN) EK(I,J)=G*(DJT(IE,1)*SN*SN+DJT(IE,2)*CS*CS) 60 CONTINUE IF(JA.EQ.2) GOTO 80 EK(3,3)=1.00 EK(3,12)=-1.00 EK(12,3)=-1.00 EK(12,12)=1.00 EK(6,6)=1.00 EK(6,9)=-1.00 EK(9

38、,6)=-1.00 EK(9,9)=1.00 80 CONTINUE RETURN END,SUBROUTINE JNT(JT,IE,NEE,IM,PA,GAM) C C Stiffness matrix of Goodmans interface element C IEJ=IE-NEE G=1.0-SS(IE)*CT(IM,4) f=a2(ie) STS(JT,IE,3)=G*G*CT(IM,5)*(A2(IE)/PA)*CT(IM,6)*GAM 555 STN(IE,3)=CT(IM,10) 566 V=SIN(BT(IEJ) W=COS(BT(IEJ) P=X(ME(IE,1),1)-

39、X(ME(IE,2),1) Q=X(ME(IE,1),2)-X(ME(IE,2),2) A=SQRT(P*P+Q*Q) DO 577 M=1,4 DO 577 N=1,4 AKA=1.0 IF(M.GT.2.AND.N.LE.2.OR.M.LE.2.AND.N.GT.2)AKA=-1.0 AKB=1.0 IF(M.EQ.N.OR.(M+N).EQ.5) AKB=2.0 G=AKA*AKB*A/6.0 I=2*M J=2*N EK(I-1,J)=G*(STS(JT,IE,3)-STN(IE,3)*V*W EK(I,J-1)=EK(I-1,J) EK(I-1,J-1)=G*(STS(JT,IE,3

40、)*W*W+STN(IE,3)*V*V) EK(I,J)=G*(STS(JT,IE,3)*V*V+STN(IE,3)*W*W) 577 CONTINUE RETURN END,三、有厚度接触面单元,Goodman单元不足： 1 单元无厚度，在受压时就会使两侧的二维单元重叠。 2 kn 任取一大值，法向位移的微小误差就会使法向应力n有较大误差。计算的n有时是不合理的,但在薄层单元本构矩阵中将法向与切向的分量分开考虑,不考虑法向和切向的耦合影响，可取 Dsn=Dns=0,Desi 薄层单元,如图接触单元长B，厚度t，在单元劲度矩阵的形成等方面与其他固体单元一样，可为4节点，8节点等参元,法向分量D

41、nn与接触面t厚度区域内材料的性质有关，也与两侧材料性质有关。可以写成,式中下标 i，g 和 st 分别表示接触区材料、岩土材料和结构材料，而1、2和3是0与1之间的系数。即Dnn是三种材料法向分量的某种加权平均值。 简单的处理是令11， 2=3=0。在动力问题中可取10.75， 2=0.25，3=0,剪切分量Dnn由接触面剪切模量构成，作接触面直剪试验，点绘与相对位移s之间的关系，剪切模量为,上式表明，单元厚度的选择对G的数值有直接的影响。研究表明宜取,B,t,土体填方工程,4,填土工程在有限元计算中与地基不同，就是在施工逐级加荷过程中，不仅荷载不断增加，而且结构本身也在逐渐扩大，即填土体在

42、加高。以下问题需要处理,1 计算网格要随施工过程而增加,可以将新填土只作为荷载作用于老的土体上，而不形成网格。在这一级末了才形成网格，参加到下一级的整体计算中去。 另一方法是新填土既作为荷载，也在该级增量中形成网格，作为整个结构的一部分来承担荷载。这两种方法均可用，以后者更普遍,2 新填土层的初始应力为0,3 新填土层完成时的应力,无法形成弹性矩阵，这时可根据单元形心到新填土层顶面的深度计算自重应力，推求弹性常数形成 D，也可直接假定弹性常数,两种计算方法： 令其等于自重应力； 新填土作为网格参加有限元计算，然后与其他老土单元一样，由节点位移求解应力。以前者更实际,传统的计算方法是假定填土一次

43、到顶，形成自下而上成三角形分布的自重应力，在自重应力作用下土层压缩产生垂直位移，填土顶面的计算压缩层最厚，计算垂直位移最大。 某一深度 z 处的的一点 A 的沉降，是h厚度的土受ABCD梯形分布的自重应力所产生的压缩量,4 填土逐级加高的计算位移，与填土完成形成荷载突然施加所产生的位移是不同的,h,z,H,A B,D C,H,对填土逐级加高的情况，如果不考虑固结等时间因素对变形的影响，即假定变形在施工中瞬时完成，则下部结构的自重不影响上部结构的变形，施工进行到什么高度，这个高度以下土重引起的位移已经发生，这个高度以上各点的垂直位移仅仅是其上土重的作用引起。 因此 A点沉降是 h 厚的土体受AB

44、ED矩形分布的应力所产生的压缩量。填土的顶面不再有荷载，也就不发生位移，所以顶面位移为0,h,z,H,A B,D E C,H,假设土体为弹性体，体积压缩系数 mv 为常量，求顶面下 z 深度处 A 点的沉降。它是 A 点以下厚度为hz的土层的压缩量。 对一次加荷来说，引起压缩的应力为全部自重应力，A点的沉降为,h,A,h,z,H,A B,D C,H,z,h,A B,D E C,h,对逐级加荷来说，当填土到A点的高度时，A点以下土体自重引起的沉降已在填土过程中完成，A点沉降为0。当填土超过 A 点高度时，A 点才有沉降，A点的沉降是A点以上土重引起的A以下土层的压缩。当填土到坝顶时，引起A点的沉

45、降的应力在压缩层hz内呈矩形分布，A点的沉降为,h,在施工逐级加荷计算整理位移成果时，对某一级荷载增量而言，荷载是一次施加的，对于该级荷载的新填土层来说，相当于对一个小的坝体作一次加荷计算，其顶面位移不为0，则各级荷载下的位移累加起来就出现阶梯状。 台阶的大小与计算分层有关，这显然不符合实际。设想对新填土层又分若干层次逐级施加，当分层无穷多时，顶面位移就是0,为了避免累计位移的台阶状，对每一新填土层可把一次荷载算出的位移，修正到分级无穷多时的位移（只对新填土层修正），近似的修正关系为,网格的处理堤坝施工,在分析开始，所有堤坝单元都不激活，在增量1分析时，施工第一层，因而在本增量开始，第一层的所

46、有单元都被激活，然后施加第一层的自重，装配总刚矩阵，进行计算。 其他各层土的填筑与此相同。通过累计每级增量分析结果可得最后计算结果。显然，结果依赖于每级填土层的厚度,layer1,layer2,layer3,1 按一组增量进行分析，在某级增量或填筑开始前，确保所有拟填筑单元未激活,2 在进行某级增量分析时，激活该级增量对应的填筑单元，并赋予适合于填筑过程中材料性质的本构模型，通常此时填筑材料的刚度较低,3 计算由于填筑材料自重引起的节点力，并加入平衡方程右边的载荷列矢量,4 组装该级增量的总刚矩阵和边界条件，解方程得该级增量位移、应变、应力,5 在下级增量开始前，改变刚填筑单元的本构模型使其能

47、表征已就位填筑材料的性质。对仅与本级填筑有关的单元节点（或新增节点）位移归零,6 进行下级增量分析,数值模拟施工过程,1 随着施工的进行，单元（新填土层）逐渐被激活，并参与计算，尚未填筑的土层单元未被激活也不参与计算（总刚与载荷矢量组装时,2 所有单元（填土层）都处于激活状态，只是在未被填筑时，其被赋予较小的刚度，但每一增量计算时所有单元都参与,数值模拟网格处理的两种方式,很多软件采用第二种处理方式,layer1,layer2,layer3,Modify displacements for new filled layers (to 770) IG=0 DO 730 IE=NF,NE DO 7

48、30 JE=1,4 I=ME(IE,JE) IF (IE.EQ.NF.AND.JE.LT.4) GOTO 720 DO 710 J=1,IG I0=IH(J) IF (I.EQ.I0) GOTO 730 710 CONTINUE 720 CONTINUE IG=IG+1 IH(IG)=I 730 CONTINUE DO 770 I00=1,IG I=IH(I00) XZ(1)=X(I,1) XZ(2)=X(I,2) CALL CYY(CY0,NC0,IC0) NC1=NC(NL) DO 740 IJ=1,NC1 740 IC1(IJ)=IC(NL,IJ,CALL CYY(CY1,NC1,IC1

49、) DZ=CY1-XZ(2) DH=CY1-CY0 IF (DH.LT.0.00001) THEN AA=1.0 GOTO 750 ENDIF AA=2.0*DZ/(DH+DZ) 750 CONTINUE DO 760 J=1,2 D1=(AA-1.0)*DU(I,J)*100.0 U(I,J)=U(I,J)+D1 UU(I,J)=UU(I,J)+D1 760 CONTINUE 770 CONTINUE,土体挖方工程,5,有限元法计算分析挖方工程，是要模拟被开挖的那部分土体去除后，地基中的应力变形状态，分析对挖坡稳定性的影响。开挖是逐步实施的，计算模拟也要分级逐步进行,一）开挖计算的主要问题,

50、1、开挖体的处理,1）保留原来的网格体系。令挖去土体的单元的材料由土体变为空气，对空气用很小的弹性模量,2）从网格中删去已挖的土体单元。这种方法更合理，但增加了计算程序中处理的麻烦，即单元数、节点数逐级减少,A B C D E F,2、应力的释放,1）直接由挖去的各层的土重求等效的节点荷载，将其反号，成为负的荷载作用在开挖面上，但这只能看作一种粗略的估算。实际上拟开挖的面上各点在开挖过程中应力状态已经发生变化,2）对开挖工程的模拟，要计算拟开挖面上的实际应力状态，将其释放为零，而成为应力自由面，这就需要根据开挖进展不断作修正计算。 当开挖到某一层时，将此时该开挖面上各点的应力分量转换为等效的节

51、点荷载，取其负值，作用在开挖面上，以抵消开挖前该面上的应力。以此作为整个地基上的荷载，进行有限元计算,第一层,第二层,A,B,由与开挖边界相邻的开挖单元的应力计算节点力,3、卸载,挖方施工时土体处于卸载状态，这里卸载是指围压降低，至于剪应力，随着开挖的进展还是增加的，仍然为加载。 在选择本构模型时要注意，模型要能反映围压降低与增加对变形的不同影响,二）开挖面上应力释放的模拟,第一级开挖模拟时，必须先求得开挖前的应力。这时一般假定土体在自重应力作用下处于 K0 状态，即平面应变条件下为,经过第一级开挖，地面不再水平，以后各级便不能再按自重应力的方法求了，而要用有限元法计算出应力。有限元法直接求得

52、的是单元应力，可用如下两种方法求边界应力,1 由单元应力内插求节点应力,2 由位移直接求节点力,在第一层开挖前，三个开挖面上的应力均为土体自重应力，三个开挖面以上土体对其下部土体所作用的等效节点力，就由自重应力求得，分别以FI，FII 和 FIII表示,第一阶段开挖，在开挖面上作用节点荷载 -FI，去除该面以上土体，作有限元计算，得出相应的地基内应力、变形。同时也可算出第二和第三开挖面上的节点力由于第一层的开挖而产生的变化量FII1和FIII1 。 FII1由紧挨在第二开挖面上方的单元劲度矩阵与第一阶段的单元节点位移增量相乘求得。此时第二开挖面上的节点力为,FII = FII + FII1,F

53、III = FIII + FIII1,同理第三开挖面上的节点力为,第二阶段开挖，在第二开挖面CEFD上作用-FII，使该面成为应力自由面。从网格中删除该面以上单元，用有限元法解得位移和应力。求第二阶段开挖在第三开挖面上所产生的节点力变化 FIII2 ，方法同前。进而求第三开挖面上当前的节点力,第三阶段开挖，就是在第三开挖面上作用节点荷载 - F”III,F”III = FIII+ FIII 2,FII1由紧挨在第二开挖面上方的单元劲度矩阵与第一阶段的单元节点位移增量相乘求得,M点上方紧挨的两个单元为3和4，分别用这两个单元的劲度矩阵与单元节点位移增量相乘，求得相应的节点力增量，其中只取M点的节

54、点力增量。将这两个单元在M点的节点力增量叠加起来，就是所求的M点的 FII1,三）水下开挖,H,A,A,H,B,G,A,H,B,G,C,F,1 去掉AH上的水压力，同时施加ABGH上的水压力,2 开挖第二层时，去掉B、G节点由于原BG面上水压力引起的节点力。因为相关部分的土体没开挖,3 BG面上B、G节点以外的其他节点力可以不管，因相关单元的开挖模拟时会间接考虑,c c 推求开挖释放等效结点荷载 c SUBROUTINE FRUE(IE,NP,W) REAL*8 RSW B=0.00 DO 10 I=1,NP 10 B=B+XZE(I,2)/FLOAT(NP) DO 40 I=1,NP IP=

55、INE(I) J=I+1 IF(J.GT.NP) J=1 JP=INE(J) I1=IWU(IP,1) I2=IWU(IP,2) J1=IWU(JP,1) J2=IWU(JP,2) IF(I1+I2).LE.0).OR.(J1+J2).LE.0) GOTO 40 DX=XZE(I,1)-XZE(J,1) DZ=XZE(I,2)-XZE(J,2) Z=0.50*(XZE(J,2)+XZE(I,2) PX=0.50*(DZ*SS(IE,1)-DX*SS(IE,3) c PY=0.50*(DZ*SS(IE,3)-DX*(SS(IE,2)+W*(Z-B) PY=0.50*(DZ*SS(IE,3)-DX

56、*SS(IE,2) IF(PX.lt.0) PX=-PX IF(PY.lt.0) PY=-PY,A=W*DX*DZ/12.00 IF(I1.NE.0) RSW(I1)=RSW(I1)-PX IF(I2.NE.0) RSW(I2)=RSW(I2)-PY IF(J1.NE.0) RSW(J1)=RSW(J1)-PX IF(J2.NE.0) RSW(J2)=RSW(J2)-PY 40 CONTINUE IF(NP.EQ.3) GOTO 60 DO 50 IC=1,4 CALL EJB(IC,2,DJ) DO 50 I=1,4 L=IWU(INE(I),2) IF(DJ.LT.0) DJ=-DJ IF

57、(L.GT.0) RSW(L)=RSW(L)-W*DJ*FN(I,IC) 50 CONTINUE GOTO 80 60 CALL EBC(1,DJ) DO 70 I=1,3 L=IWU(INE(I),2) 70 IF(L.GT.0) RSW(L)=RSW(L)-W*DJ/3.0 80 CONTINUE RETURN END c,空气单元,平面问题的地下结构单元,6,1. 墙,混凝土墙： 实体单元,钢板桩、混凝土桩或其他钢与板组合墙：梁单元,2. 地锚,实体单元、 弹簧单元、 接触面单元,弹簧单元,弹簧单元,A,B,C,膜单元,膜单元,A,B,C,弹簧单元,A,B,C,弹簧单元、膜单元,膜单元,

58、3. 墙与支撑的连接,在平面应变分析中，挡土墙与支撑的连接有以下方式,Shear =0 Moment = 0,Tied A and B in Horizontal direction,Axial thrust 0,Solid elements,Tied A and B in Horizontal direction,Beam/bar elements,A B,简单连接,Moment = 0,Axial thrust 0 shear 0,Tied A and B inHorizontal and Vertical direction,Solid elements,Solid elements,B

59、eam/bar elements,A B,销连接,Axial thrust 0 shear 0 moment 0,Solid elements,Beam/bar elements,完全连接,4. 耦合分析中的结构单元,在耦合固结分析中，挡土墙如图a。用梁单元模拟墙如图b，由于实体单元的节点为普通节点，在固结分析中实际意味墙是透水的，水可自由流过。这与实际不符,如果要使墙为不透水，可以在墙的一边设置接触面单元如图c。如果这些接触面为非固结单元，他们将在墙的某一边提供一个不透水的隔断,a b c,浅基础与桩基础的有关计算,7,面基础,浅基础,1 柔性基础,在传统土力学中，分析面基础时，通常假设为两

60、种极端情况，即理想柔性基础和理想刚性基础。在数值计算时如果也用这种假设，则在计算过程中可以不包括基础本身，而以边界条件形式考虑,Y(v,X(u,载荷控制,位移控制,不可能,光滑,粗糙,Fy应用 Fx=0 （一般不设定边界力时，软件默认为0,Fy 应用 u = 0,2 刚性基础,Y(v,X(u,光滑,粗糙,Fy应用 Fx=0 v 捆绑,Fy应用 u = 0 v 捆绑,载荷控制,Y(v,X(u,光滑,粗糙,v 应用 Fx=0,v 应用 u = 0,位移控制,基础荷载,粗糙,接触面单元,位移控制,v 应用,3 浅基础,接触面单元,4 基础掏土纠偏,目前塔北侧沉降1m余，南侧下沉近3m，沉降差达1.8