【2019年整理】岩土工程数值计算.ppt

1、岩土工程数值计算二零一一年三月参考文献参考文献计算土力学 上海科学技术出版社 朱百里土工计算机分析 中国建筑工业出版社 龚晓南土工数值计算 中国铁道出版社 钱家欢土工原理与计算 水利电力出版社 钱家欢岩土工程有限元分析理论与应用 科学出版社 谢康和Finite element analysis in geotechnical engineering David M.Potts 有限元法的基本知识有限元法的基本知识 0.1 1 自由度的编号自由度的编号 N=n+1 READ(1,*)(IB3(I),I=1,NP3)READ(1,*)(LH(I),I=1,NP5)DO 60 I=1,NP3 IS=

2、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 100 IP3(I,J)=IP3(I,J)+N3 N3=N3+1100 CONTINUE80 CONTINUE213 0 1节点号 X 方向约束信息Y 方向约束信息121 0 节点号 0 代表孔压2 2 半带宽 半带宽 dmdm dm=b

3、+1 b为相邻自由度编号的最大差值（对应节点在同一单元）nnmmkkkkkkkkkkkkkkkkkK99887766555453524544353325221100dm 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(I

4、CNP)=ICNP-ICNQ20 CONTINUEHALF BAND WIDTH-1 RELATED TO Ith FREE DEGREE3 3 总刚的组装总刚的组装 节点平衡法12345每一单元节点力与节点位移之间的关系为 22222RFFFFe321333231232221131211321kkkkkkkkkFFF5214212213212211211kkkkkkF3212211213212211212kkkkkkF2213215213212211213kkkkkkF根据外荷载与节点力平衡条件 直接刚度法先把每个单元的单刚阶数扩大成总刚阶数，把单刚中按局部编码的子块搬到总刚中相应的总体编码

5、的位置中去，余下部分用零子块填充。123456 32100000000000000000000000000654321333231232221131211kkkkkkkkkK1 2 3 4 5 61 2 3 4 4 总刚的存储总刚的存储 一维变带宽存储由于整体刚度矩阵具有对称性、稀疏性和带状性，采用变带宽下三角一维存储。887877686766565546454435343324232214131211kkkkkkkkkkkkkkkkkkkkk887877686766565546454435343324232214131211kkkkkkkkkkkkkkkkkkkkk DO 20 I=1,NP

6、 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将某一单元的 Ke 集整到劲度矩阵 K中 SKYLINE(轮廓线法)按列存储刚度矩阵上三角区必要部分，按行存储下三角区中必要部分。对刚度矩阵中出现少数非常长的列的情况下，存储要求不会剧烈增加，

7、很容易利用向量点积例行程序。8878776867665856554846454438353433282423221814131211kkkkkkkkkkkkkkkkkkkkkkkkkkK11 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+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

8、)-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 CONTINUE 345 CONTINUE ENDIF SL(NR)=SL(NR)+ZY(NRCC)350 CONTINUE 355 CONTINUE(NCN,NR)(NR,NCN)5 5 边界条件的引边界条件的引入入 划 0 置 1

9、 法：处理 ui=0 约束 乘大数法：处理 ui=R 约束将总刚相应的主对角元素改为 1，将对应的行、列其它元素改为 0。将荷载向量中相应的元素改为 0。或预先将每个节点的方程编号，已知位移的节点不编号。将总刚相应的主对角元素置一大数，将荷载向量中相应的元素改为该大数乘 R。初等变换法：处理 两个变量有确定关系如轴对称问题，可以取一夹角为扇形区计算。在斜对称 AB边上的点只能沿 AB上移动，所以点 C上的点的位移有使用初等变换法消除一个相关方程，若消除对应的方程，可以将 所在的行乘 tg 加到 u 所在的行，并将 所在的列也乘 tg 加到 u 所在列上去。tguccJP(1):NODE NUM

10、BERJP(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 选择单元类型、本构模型，并输入本构参数5 计算结果提取、分析等岩土工程问题分析方法岩土工程问题分析方法 1.1.1 1.1 岩土工程问题控制方程的建立岩土工程问题控制方程的建

11、立0uzyxzxyxx 0vzyxzyyxy 0gwzyxzyzxz 1 1 土体平衡方程土体平衡方程2 2 土体本构方程土体本构方程 Ddzzppdxxppgzxwwzpp3 3 土体几何方程土体几何方程xuxyvyzwzxvyuxyywzvyzzuxwzx4 4 土体有效应力原理土体有效应力原理pxxpyypzz5 5 孔隙流体平衡方程孔隙流体平衡方程0ukvgxpwxxw 0vkvgypwyyw 0gwkvgzpwwzzw 6 6 渗流连续方程渗流连续方程zwyvxutzvyvxvzyxdzzzzdzxxxgzxw7 7 总控制方程总控制方程uxpzxwddyxvddzudyudxud

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

13、形和稳定性分析。0)()(2551324412225522442211zxwddyxvddzudyudxud0)()(2551324412225522112244zywddyxuddzvdyvdxvd0)()(2551325513223322552255zyvddyzuddzwdywdxwd2 2 有效应力分析法及其控制方程有效应力分析法及其控制方程在有效应力分析法中，土体的有效应力和孔压被严格区分，并将土骨架变形与孔隙水的渗透同步考虑。因此，有效应力分析法较能更真实地反映土体的自身特性，能更合理地计算土体对载荷的响应，应用范围更广。有效应力分析法尚需要有效应力原理和连续性方程。总应力法中只有

14、位移变量且仅与空间有关，而有效应力法中还有孔压变量，而且与空间和时间均有关。0)()(2551324412225522442211xpzxwddyxvddzudyudxud0)()(2551324412225522112244ypzywddyxuddzvdyvdxvd0)()(2551325513223322552255gzpzyvddyzuddzwdywdxwdzwyvxutggzpkypkxpkwwzyx2222223 3 总应力分析法总应力分析法和和有效应力分析法有效应力分析法关系关系总应力法是有效应力法中当孔压 p=0时的特殊形式。在有效应力分析中如果采用与总应力分析相同的土工参数，并

15、令孔压 p=0，所得结果即为总应力分析结果。总应力分析一般采用土体得不排水指标，由此进行的是加荷瞬时或短期应力和变形分析。但也可采用土体的排水指标，此时进行的是最终或长期应力和变形分析。当进行线弹性分析并采用排水指标，总应力分析得到的结果为有效应力分析的最终（孔压消散完毕、主固结完成）结果。4 4 不排水孔压计算不排水孔压计算 f 000ffffpppDffD fDDD333300IIKDef有效应力原理其中fsffpKnpKn1vv)(ezyxefKKpsfeKnKnK11孔压增量不仅使孔隙流体压缩，也会引起土颗粒的体积压缩。相应的有效应力增量也会引起土颗粒的体积变化，然而，由于有效应力必须

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