岩土工程数值计算2011

岩土工程数值计算 二零一一年三月 参考文献 计算土力学上海科学技术出版社朱百里 土工计算机分析中国建筑工业出版社龚晓南 土工数值计算中国铁道出版社钱家欢 土工原理与计算水利电力出版社钱家欢 岩土工程有限元分析理论与应用科学出版社谢康和 FiniteelementanalysisingeotechnicalengineeringDavidM Potts 有限元法的基本知识 0 1自由度的编号 N n 1 READ 1 IB3 I I 1 NP3 READ 1 LH I I 1 NP5 DO60I 1 NP3IS IB3 I CALLHDIR IS DO60J 2 3IP3 JP 1 J 1 JP J 60CONTINUEDO70I 1 NP5IS LH I 10II LH I IS 10IP3 IS 3 II70CONTINUEN3 0DO80I 1 NPDO100J 1 3IF IP3 I J EQ 0 GOTO100IP3 I J IP3 I J N3N3 N3 1100CONTINUE80CONTINUE 21301 节点号X方向约束信息Y方向约束信息 1210 节点号0代表孔压 2半带宽dm dm b 1 b为相邻自由度编号的最大差值 对应节点在同一单元 dm IJK 3NN 12DO10I 1 4IVP IVEN I DO10J 1 IJKIJ IJK I 1 JICN IJ IP IVP J 10CONTINUEDO20J 1 NNICNP ICN J DO20K 1 NNICNQ ICN K IF ICNP EQ 0 OR ICNQ EQ 0 GOTO20IF IA3 ICNP GT ICNP ICNQ GOTO20IA3 ICNP ICNP ICNQ20CONTINUE HALFBANDWIDTH 1RELATEDTOIthFREEDEGREE 3总刚的组装 节点平衡法 1 2 3 4 5 每一单元节点力与节点位移之间的关系为 根据外荷载与节点力平衡条件 直接刚度法 先把每个单元的单刚阶数扩大成总刚阶数 把单刚中按局部编码的子块搬到总刚中相应的总体编码的位置中去 余下部分用零子块填充 1 2 3 4 5 6 123456 123 4总刚的存储 一维变带宽存储 由于整体刚度矩阵具有对称性 稀疏性和带状性 采用变带宽下三角一维存储 DO20I 1 NPDO20J 1 JANX JA I 1 JM INE I 20MA NX IWU M J M NP JADO60I 1 MK MA I IF K LE 0 GOTO60DO40J 1 ML MA J IF L EQ 0 OR L GT K GOTO40II IDK K L KTK II TK II EK I J 40CONTINUE60CONTINUE 将某一单元的 K e集整到劲度矩阵 K 中 SKYLINE 轮廓线法 按列存储刚度矩阵上三角区必要部分 按行存储下三角区中必要部分 对刚度矩阵中出现少数非常长的列的情况下 存储要求不会剧烈增加 很容易利用向量点积例行程序 K11k12k22k13k23k33k14k24k34 DO355K 1 4NRCC 3 K 1 NR NQ NOD K IX 1DO350M 1 2 IFLOWNRCC NRCC 1NR NR 1IF ICODE LT 2 THENLENNC LOCC NR LOCC NR 1 1 2DO345L 1 4NCCC 3 L 1 NCN NQ NOD L IX 1DO344N 1 2 IFLOWNCCC NCCC 1NCN NCN 1IF NR LT NCN NN LOCC NCN NCN NRIF NR EQ NCN NN LOCC NR IF NR GT NCN NN LOCC NR LENNC NR NCNS NN ISHIFT S NN ISHIFT C1 NRCC NCCC 344CONTINUE345CONTINUEENDIFSL NR SL NR ZY NRCC 350CONTINUE355CONTINUE NCN NR NR NCN 5边界条件的引入 划0置1法 处理ui 0约束 乘大数法 处理ui R约束 将总刚相应的主对角元素改为1 将对应的行 列其它元素改为0 将荷载向量中相应的元素改为0 或预先将每个节点的方程编号 已知位移的节点不编号 将总刚相应的主对角元素置一大数 将荷载向量中相应的元素改为该大数乘R 初等变换法 处理两个变量有确定关系 如轴对称问题 可以取一夹角为扇形区计算 在斜对称AB边上的点只能沿AB上移动 所以点C上的点的位移有 使用初等变换法消除一个相关方程 若消除 对应的方程 可以将 所在的行乘tg 加到u所在的行 并将 所在的列也乘tg 加到u所在列上去 JP 1 NODENUMBERJP 2 X DIRECTIONRESTRAININFORMATIONJP 3 Y DIRECTIONRESTRAININFORMATION JP 1 I3 100K3 I3 JP 1 100JP 2 K3 10K3 K3 JP 2 10JP 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总应力分析法及其控制方程 总应力分析法与一般固体力学相同 从应用上讲 一般用于不考虑渗流固结的情况 如饱和粗粒土地基 透水性土料组成的土坝路堤的应力和变形分析以及饱和软粘土地基短期变形和稳定性分析 2有效应力分析法及其控制方程 在有效应力分析法中 土体的有效应力和孔压被严格区分 并将土骨架变形与孔隙水的渗透同步考虑 因此 有效应力分析法较能更真实地反映土体的自身特性 能更合理地计算土体对载荷的响应 应用范围更广 有效应力分析法尚需要有效应力原理和连续性方程 总应力法中只有位移变量且仅与空间有关 而有效应力法中还有孔压变量 而且与空间和时间均有关 3总应力分析法和有效应力分析法关系 总应力法是有效应力法中当孔压p 0时的特殊形式 在有效应力分析中如果采用与总应力分析相同的土工参数 并令孔压p 0 所得结果即为总应力分析结果 总应力分析一般采用土体得不排水指标 由此进行的是加荷瞬时或短期应力和变形分析 但也可采用土体的排水指标 此时进行的是最终或长期应力和变形分析 当进行线弹性分析并采用排水指标 总应力分析得到的结果为有效应力分析的最终 孔压消散完毕 主固结完成 结果 4不排水孔压计算 有效应力原理 其中 孔压增量不仅使孔隙流体压缩 也会引起土颗粒的体积压缩 相应的有效应力增量也会引起土颗粒的体积变化 然而 由于有效应力必须通过颗粒接触 但接触面积很小 导致体积变化也很小 如果忽略这种体积变化 则总的体积变化为 对于饱和土体 Ks Kf都比土骨架模量大很多 在不考虑具体值时 此时准确值已不重要 可假设Kf Ks 得 由于Ks比土骨架模量大很多 如孔隙流体压缩性较大 以致Ks kf 则 对于排水分析 取Ke 0 加载过程中孔压不变 对于各向同性线弹性土体 进行不排水分析时 必须设置Ke 据经验 对于饱和土体 只要ke足够大 土体对ke的实际大小并不敏感 但Ke太大时 可能导致数值不稳定 即不排水泊松比接近0 5 专家建议设置ke Kskel 在100和1000之间 Kskel是土骨架的体积模量 1 3岩土工程问题的边界条件 1固结分析中的边界条件 对节点位移和孔压已知的情况 可以有两种处理方法 1 仍给以自由度编号 在解方程时把解得的位移值用已知值置换 2 已知值对应的自由度均编为0 把代数方程区分为 式中 为已知的边界值 计算中只对x1建立方程 即 Displacement porepressureDOF DisplacementDOF 其中 A12 x2 项在对各单元进行计算时即可求出 在多数情况下 已知的边界值为0 上式简化为 2渗流问题的边界条件 岩土工程种进行渗流计算的主要目的是求得渗流量和孔隙水头或水力坡降 渗流边界可分为两类 即已知水头边界和已知流量边界 分别可以表示如下 在渗流计算中 对于自由面渗流问题 其边界条件需特殊处理 1 有压渗流 有压渗流的计算域边界和边界条件是固定的 可以很方便地应用有限元法计算 如图 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 自由面渗流 在堤坝 渠道和地下水渗流计算中 自由面或浸润面事先是未知的 因此计算域本身是未知的 早期的研究工作是用试算法不断改变自由面位置 直到计算结果满足如下自由面上的边界条件为止 1 2 然后用新坐标进行下一轮计算 直到自由面上所有节点满足h z 精度要求 可取最大水头差的1 计算中把自由面当作不透水面 所以式 2 自动满足 而式 1 一般不满足 如果某一节点在第i次计算中假定的节点高程为zi 求得的节点水头为hi 不满足式 1 则利用下式计算新的坐标 具体计算时把渗流域分成非饱和的负压区 饱和的正压区和过渡区三部分 并把渗透系数随孔压的变化简化成三段折线 假定负压区的渗透系数很小 等于饱和区的1 过渡区内渗透系数随孔压p g h z 直线变化 此时下图边界条件为 h1 h0 gf a e d bc ab和ag bc和fg cd de和ef 另一种方法是把自由面以上的非饱和区和自由面以下的饱和区同时计算 计算所得的零孔隙压力线即为满足h z条件的自由面 3 不稳定渗流 随时间变化的渗流称为不稳定渗流 随动力荷载迅速变化的渗流问题需结合土骨架的变形同时计算 水库中水位降落时的渗透力将对土坝上游边坡造成不利影响 是不稳定渗流计算的主要对象 不考虑水体的压缩性 不稳定渗流控制方程仍如前所述 只是边界条件随时间变化 特别在自由面渗流中 自由面也是变化的 此时自由面除应满足h z外 还应满足下列边界条件 式中ne为有效孔隙率 相当于单位土体孔隙中自由水的含量 许多文献将上式表示为 此式的物理含义是 进入自由面的水量恰为自由面变动区的土体 图中阴影区 内所含的水量 ne为给水度 把写成差分形式 h0为上一时段的水头值 并代入有限元插值公式 则某一节点i的渗出量为 nf为自由面节点数 如此 前述有限元公式可以改写为 如此 前述有限元公式可以改写为 nf2 nf nf1 nf1为自由面上水头已知节点数 nt中仍包括自由面节点nf2 nb1为包括nf1在内的所有水头已知节点 由于自由面上还需满足h z的条件 所以对每一时段 t还需进行自由面位置的迭代计算 f 真空堆载预压分析 2 2 1施工工艺 先打袋装砂井 直径约7cm 间距约1m 上面铺50cm的砂垫层 其中埋有真空滤管及与其连接的主管 然后在加固区四周开挖深约1m的沟槽 铺上塑料薄膜 将其周边埋入沟中 并将主管伸出薄膜与真空设备相连 密封沟 抽真空系统 主管 支滤管 1234 1 2 3 4 2 2控制方程 在平面上用有限单元法离散化 时间上用差分法分段以后 BIOT固结理论可以表示为 其中 和为时段内j节点的水平 垂直位移和孔压增量 为 i结点的水平 垂直向荷载增量 pj0为上一时段末j结点的孔隙压力 为差分常数 此处用2 3 等为方程式的系数 反映j结点对i结点的影响 其中 和 与选用的单元形状有关 除单元形状外 还与土的渗透系数有关 而 和除单元形状外和 还与土骨架的变形性质有关 作为一个边值问题 土体内部发生固结变形过程是边界条件改变的结果 边界条件可以分为荷载边界 位移边界和水头边界 水量边界几种 真空预压与加荷预压只是在边界条件上有所不同 而控制方程式与解题方法是完全一样的 在加荷压缩问题中 只有边界荷载 应力控制 或边界位移 应变控制 发生变化 水流边界条件不变 在真空抽水问题中 则只有水流边界条件发生变化 即真空作用点的水头下降 从而形成由内向外的水力坡降 显然 两者同时变化的问题 真空联合堆载 也可同样求解 说明 一般固体力学问题的属性通常仅指应变或应力属性 而固结问题的属性则包括应变和渗流属性 在固体力学中 平面问题指应变或应力是二维的问题 但对固结而言 当变形和渗流均为二维时才属于平面问题 2 3固结问题属性 在一般的实际荷载作用下 天然软土地基或人工处理后地基固结问题通常属于变形和渗流均为三维的空间问题 但在长条形荷载作用下 天然软土地基固结问题属于平面变形 平面渗流固结问题 即平面问题 而经人工处理后的地基固结则往往属于平面变形 空间渗流固结问题 为了编程和计算方便 常将砂井固结问题作为一般的平面问题 平面变形 平面渗流 分析 将砂井简化为沿荷载长度方向连续 垂直于计算断面 的砂墙 并通过改变砂井和土体的渗透系数等来使原砂井地基与简化后的砂墙地基等效 2 4砂井 排水板 的处理 第 种方法 将竖向排水通道视为渗透性很强的介质 由于砂井的间距很小 而为了得到合理的孔隙压力分布 两井之间至少应布置三排结点这样势必使结点数成倍增加 增加计算工作量 因此还需要放大砂井的间距 但这两种变换应保证变换前后主要基本量 如固结度或同一深度处的平均孔压 保持不变的前提下进行 这可以通过在平面应变有限元计算时调整砂井地基的实际渗透系数来实现 把土的水平向渗透系数按井距放大倍数的平方放大 以保持水平向相对的渗径长度不变 这样的处理方法虽然比较粗糙 但经验表明计算结果大体上还是能反映实际情况的 根据砂墙地基双向应变 双向渗流等应变固结理论解与巴隆轴对称固结理论解的比较 在固结度或平均孔压不变的条件下 砂井地基与砂墙地基的等效可通过调整渗透系数得到 其等效计算公式为 Kxp Dx kra Kzp Dz kza Dx Dz为水平向和竖向渗透系数的调整系数 kxp kzp分别为砂墙地基的水平和垂直渗透系数 kra kza分别为砂井地基的水平和垂直渗透系数 实际计算时 用塑料排水板代替砂井 塑料排水板的当量直径采用下式计算 式中 Dp为塑料排水板当量直径 cm 为换算系数 无实验资料时可取0 75 1 00 b为塑料排水板宽度 为塑料排水板厚度 塑料排水板有效排水区直径按下式计算 式中 de为塑料排水板有效区直径 S为塑料排水板间距 de 1 05S 排水板与砂井的等效 在进行有限元计算时 根据实测膜下真空度值 将密封膜覆盖的地面边界单元结点的孔压设为负的稳定孔压 膜下真空度 在砂墙的位置 沿纵向与土接触的边界单元结点的孔压根据实测真空度调节为稳定负的孔压 从地面向下由大到小按线形变化设置 砂墙位置的底部负孔压为零 这与观测的孔隙水压力基本相符合 砂墙位置土的渗透系数及各种参数则跟周围土层相同 第 种方法 将竖向排水通道视为负压边界条件 岑仰润 真空预压加固地基的试验及理论研究 D 浙江大学 2003 设砂垫层和砂井中所有结点的孔隙水应力为常量 砂垫层以外地表处的孔隙水应力为零 其他边界的孔隙水应力皆未知 设地表面自由变形 S1 透水边界S2 负压边界S3 自由边界S4 不透水边界S5 对称边界 由于竖向排水体井阻与涂抹作用的存在 抽真空作用在竖向排水体中形成的负压分布是复杂的 一些学者提出的真空预压加固地基时竖向排水体中的负压分布模式相差比较大 根据真空 堆载联合预压法的加固机理可知 真空荷载主要是作用在地表 并经由塑料排水板向下传递 使土体中的孔压与地表及排水板中的孔压存在压差而渗出 从而发生固结 由于测量手段的局限 塑料排水板中的孔压变化规律尚不完全清楚 但是由于井阻等因素 塑料排水板中孔压的增长没有地表真空度那么迅速和明显 由于塑料排水板中的孔压变化是由地表真空荷载变化引起 因此 本文将地表的真空荷载作为已知边界条件 而将砂井中的孔压作为未知量处理 即将加固区地表的孔压变化值作为边界条件 堆载荷载作为一般荷载施加在加固区地表 第III种方法 砂井孔压未知 三维 孔压边界取为 砂垫层中所有结点的孔隙水压力为负的真空压力 80kPa 砂垫层以外的地基表面孔隙水压力为0 其它边界的孔压未知 这是因为真空预压时 首先降低密封膜下砂垫层中的孔隙水压力 形成膜下真空度 并通过排水通道向下传递 使得土体内部各点与排水通道及砂垫层中形成压差 从而 发生由土中向边界渗流 而在固结过程中 竖向排水通道的孔压是随时间变化的 孔隙水压力在加固过程中保持不变的只有砂垫层和砂垫层以外的地基表面 计算真空预压时 让加固区表面各点的孔压按线性从0减少到 80kPa 然后保持不变 在堆载时 将填土折算成等效结点荷载 按实际加载曲线施加上去 真空卸载时 让加固区表面各点的孔压按线性从 80kPa增加到0 竖向边界面 水平位移约束 竖向位移无约束 不透水边界 基底 竖向 水平向位移约束 不透水边界 地面 加固区和非加固区位移自由 非加固区孔压为0 加固区表面在抽真空时孔压大小为真空度 没有抽真空时该处的孔压也为0 真空度的变化和堆载的加载曲线 如图所示 具体而言 其约束条件为 假定在进行真空预压前 各单元结点孔隙水应力皆为大气压力 为解方程的方便 我们以下都取相对压力 即设pa 0 于是有 式中 0 非边界结点的孔隙水应力初始值 设在抽气的瞬时 砂井及砂垫层中立即达到相同的负压值 且保持不变 假定结点初始位移为零 即 0 0 式中 0 初始结点位移 0 0 2 5初始条件 1真空荷载处理 在砂垫层中抽真空时 负压作用在孔隙水上 使土中水部分流出土体 从而孔隙水应力降低 土骨架应力增加 而结点荷载始终为零 所以荷载向量 R 0 在解有限单元公式的方程组时 将负压边界 砂垫层 砂井 上的结点孔隙水应力的已知值乘以总系数矩阵中相应系数 移到方程右端作为荷载向量 2 6荷载处理 2荷堆处理 对填土荷载的处理一般有2种方法 1 把填土当成新增单元参加计算 2 将填土折算成等效结点荷载进行计算 前者能够模拟填土与地基的相互作用 计算填土的应力和应变 分析粘性填土的变形和开裂问题具有明显的优势 但不能模拟线性加载过程 而后者可以方便地模拟实际的加载过程 本文中采用第二种方式 填土荷载按施工设计分10级施加 计算前将加载区域都设置为null 空单元 随后分别按时间顺序给所需加载的土层赋值 从而实现分级加载 计算模型也相应地分10步完成 真空 堆载联合预压加固路堤FLAC数值分析 真空预压加固地基 均在地基中设置竖向排水体 竖向排水体的存在影响真空预压加固地基的效果和地基的固结速度 竖向排水体一般采用袋装砂井和塑料排水板 早期多采用袋装砂井 如今基本上均采用塑料排水板 当采用袋装砂井时 赋予砂井单元砂料的渗透系数 即己经考虑了井阻作用 若要再考虑涂抹作用 则将砂井单元周围土体划分为涂抹单元 涂抹区土体的强度和渗透系数适当折减以考虑涂抹作用 当采用塑料排水板时 可以将塑料排水板参照上述方法进行处理 从实际分析情况看 这样往往会过高估计负压沿竖向排水体的传递能力 2 7竖向排水体的界定 考虑到竖向排水体井阻与涂抹作用的存在 抽真空作用在竖向排水体中形成一定的负压分布 也可认为地基是在竖向排水体中负压作用下固结的 因此也可以预先设定竖向排水体中负压分布 从实际应用看 正确预先估计竖向排水体中负压分布模式是分析的关键 竖向排水体的界定 参考文献 2 软土地基真空排水预压的固结变形分析 岩土工程学报 1986 NO 3 1 用真空加固软土地基的机制与计算方法 岩土工程学报 1986 NO 2 3 刘汉龙 真空 堆载预压处理高速公路软基的有限元计算 岩土力学 2003 No 6 4 岑仰润 真空预压加固地基的试验及理论研究 D 浙江大学 2003 5 真空 堆载联合预压加固路堤FLAC数值分析 6 雷鸣 真空预压数值分析中竖向排水通道的讨论 铁道学报 2009 土与结构接触面单元 3 一 概述 对于任何土和结构相互作用情况 结构与土之间可能会发生相对运动 使用接触面单元 能模拟土与结构的接触边界条件 其优点是可以改变接触面的本构特性 允许土与结构之间的相对运动即滑动和分离 已有很多方法用于模拟土与结构接触面的非连续特性 当使用连续单元时 由于变形协调 会禁止土与结构接触面的相对运动如图 有限元法的节点协调条件将约束相邻的结构与土单元使其一起运动 这显然与实际不符 薄的连续单元 标准的本构关系 PandeandSharm1979 Griffiths1985 2 联动单元 仅考虑相对节点之间的连接 通常采用独立的弹簧连接Hermann1978 Franketal1982 3 接触面或节理单元 零厚度或有限厚度 Goodmanetal1968 Ghaboussietal1973 CarolandAlonso1983 Wilson1977 Desaietal1984 Beer1985 4 混合法 上述接触面的处理方法中 以零厚度接触面单元 Goodman单元 应用最为普遍 后又被帝国理工学院Day1990予以发展 土和结构分开模拟 通过约束方程联系 使其在接触面上保持位移与力的协调性 FrancavillaandZienkiewicz1975 Katon1983 Laiandbooker1989 二 无厚度接触面单元 Goodman等人提出一种接触面单元如图 由两片长度为L的接触面1 2和3 4组成 两片接触面间假想为无数微小的法向的和切向的弹簧连接 在受力前两接触面完全吻合 即单元没有厚度只有长度 是一维单元 一个接触面单元共4个节点如图 12 43 顶 底 z x 12 43 顶 底 z x 取线性位移模式 可将接触面上任一点的位移用节点位移表示 u3 v3 接触面单元内各点的相对位移为 式中 令 由虚位移原理可得 劲度系数 应力与相对位移的关系 取值 对于法向劲度系数 当接触面受压时 为了模拟两边二维单元不会在接触面处重叠 应取一很大的数值 如kn 108kN m3 可使相互嵌入的相对位移小到可略去 当接触面受拉时 因为接触面上不能承受拉应力 则令kn为很小的值 如取kn 10kN m3 以使算出的拉应力可忽略 拉开 重叠 对于切向劲度系数 可由直剪试验确定 CloughandDuncan摩擦试验表明 剪应力与相对剪切位移近似于双曲线关系 相应的表达式为 CloughandDuncan把初始剪切劲度和极限剪应力值与 n的关系表示为 切线剪切劲度系数为 ksi u 考虑孔压的单元刚度矩阵 K3 3 1K6 6 1K9 9 1K12 12 1 K3 12 1K6 9 1K9 6 1K12 3 1 推求接触面单元的劲度矩阵 K ecSUBROUTINEFJK IE JA COMMON PMEE PME 19 XZEE XZE 4 2 DDJJ DJT 80 10 COMMON EEKK EK 12 12 ER 12 DO40I 1 12DO40J 1 1240EK I J 0 00PA 98 10G 1 00 DJT IE 4 PME 4 F DJT IE 10 IF F LT 0 50 PA F 0 50 PADJT IE 1 G G PME 5 F PA PME 6 0 10 PAIF DJT IE 4 GT 0 99990 DJT IE 1 0 10 PADJT IE 2 PME 7 0 10 PAIF DJT IE 10 LE 0 1E 4 PA DJT IE 2 0 10 PASN 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 DO60M 1 4DO60N 1 4AKA 1 00 IF M GT 2 AND N LE 2 OR M LE 2 AND N GT 2 AKA 1 00AKB 1 00IF M EQ N OR M N EQ 5 AKB 2 00G AKA AKB A 6 00I JA M 2 JAJ JA N 2 JAEK I 1 J G DJT IE 1 DJT IE 2 SN CSEK 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 60CONTINUEIF JA EQ 2 GOTO80EK 3 3 1 00EK 3 12 1 00EK 12 3 1 00EK 12 12 1 00EK 6 6 1 00EK 6 9 1 00EK 9 6 1 00EK 9 9 1 0080CONTINUERETURNEND SUBROUTINEJNT JT IE NEE IM PA GAM CCStiffnessmatrixofGoodman sinterfaceelementCIEJ IE NEEG 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 GAM555STN IE 3 CT IM 10 566V SIN BT IEJ W COS BT IEJ P X ME IE 1 1 X ME IE 2 1 Q X ME IE 1 2 X ME IE 2 2 A SQRT P P Q Q DO577M 1 4DO577N 1 4AKA 1 0IF M GT 2 AND N LE 2 OR M LE 2 AND N GT 2 AKA 1 0AKB 1 0IF M EQ N OR M N EQ 5 AKB 2 0G AKA AKB A 6 0I 2 MJ 2 NEK I 1 J G STS JT IE 3 STN IE 3 V WEK I J 1 EK I 1 J EK I 1 J 1 G STS JT IE 3 W W STN IE 3 V V EK I J G STS JT IE 3 V V STN IE 3 W W 577CONTINUERETURNEND 三 有厚度接触面单元 Goodman单元不足 1单元无厚度 在受压时就会使两侧的二维单元重叠 2kn任取一大值 法向位移的微小误差就会使法向应力 n有较大误差 计算的 n有时是不合理的 但在薄层单元本构矩阵中将法向与切向的分量分开考虑 不考虑法向和切向的耦合影响 可取Dsn Dns 0 Desi薄层单元 如图接触单元长B 厚度t 在单元劲度矩阵的形成等方面与其他固体单元一样 可为4节点 8节点等参元 法向分量 Dnn 与接触面t厚度区域内材料的性质有关 也与两侧材料性质有关 可以写成 式中下标i g和st分别表示接触区材料 岩土材料和结构材料 而 1 2和 3是0与1之间的系数 即 Dnn 是三种材料法向分量的某种加权平均值 简单的处理是令 1 1 2 3 0 在动力问题中可取 1 0 75 2 0 25 3 0 剪切分量 Dnn 由接触面剪切模量构成 作接触面直剪试验 点绘 与相对位移 s之间的关系 剪切模量为 上式表明 单元厚度的选择对G的数值有直接的影响 研究表明宜取 B t 土体填方工程 4 填土工程在有限元计算中与地基不同 就是在施工逐级加荷过程中 不仅荷载不断增加 而且结构本身也在逐渐扩大 即填土体在加高 以下问题需要处理 1计算网格要随施工过程而增加 可以将新填土只作为荷载作用于老的土体上 而不形成网格 在这一级末了才形成网格 参加到下一级的整体计算中去 另一方法是新填土既作为荷载 也在该级增量中形成网格 作为整个结构的一部分来承担荷载 这两种方法均可用 以后者更普遍 2新填土层的初始应力为0 3新填土层完成时的应力 无法形成弹性矩阵 这时可根据单元形心到新填土层顶面的深度计算自重应力 推求弹性常数形成 D 也可直接假定弹性常数 两种计算方法 令其等于自重应力 新填土作为网格参加有限元计算 然后与其他老土单元一样 由节点位移求解应力 以前者更实际 传统的计算方法是假定填土一次到顶 形成自下而上成三角形分布的自重应力 在自重应力作用下土层压缩产生垂直位移 填土顶面的计算压缩层最厚 计算垂直位移最大 某一深度z处的的一点A的沉降 是h厚度的土受ABCD梯形分布的自重应力所产生的压缩量 4填土逐级加高的计算位移 与填土完成形成荷载突然施加所产生的位移是不同的 h z H AB DC H 对填土逐级加高的情况 如果不考虑固结等时间因素对变形的影响 即假定变形在施工中瞬时完成 则下部结构的自重不影响上部结构的变形 施工进行到什么高度 这个高度以下土重引起的位移已经发生 这个高度以上各点的垂直位移仅仅是其上土重的作用引起 因此A点沉降是h厚的土体受ABED矩形分布的应力所产生的压缩量 填土的顶面不再有荷载 也就不发生位移 所以顶面位移为0 h z H AB DEC H 假设土体为弹性体 体积压缩系数mv为常量 求顶面下z深度处A点的沉降 它是A点以下厚度为h z的土层的压缩量 对一次加荷来说 引起压缩的应力为全部自重应力 A点的沉降为 h A h z H AB DC H z h AB DEC h 对逐级加荷来说 当填土到A点的高度时 A点以下土体自重引起的沉降已在填土过程中完成 A点沉降为0 当填土超过A点高度时 A点才有沉降 A点的沉降是A点以上土重引起的A以下土层的压缩 当填土到坝顶时 引起A点的沉降的应力在压缩层h z内呈矩形分布 A点的沉降为 h 在施工逐级加荷计算整理位移成果时 对某一级荷载增量而言 荷载是一次施加的 对于该级荷载的新填土层来说 相当于对一个小的坝体作一次加荷计算 其顶面位移不为0 则各级荷载下的位移累加起来就出现阶梯状 台阶的大小与计算分层有关 这显然不符合实际 设想对新填土层又分若干层次逐级施加 当分层无穷多时 顶面位移就是0 为了避免累计位移的台阶状 对每一新填土层可把一次荷载算出的位移 修正到分级无穷多时的位移 只对新填土层修正 近似的修正关系为 网格的处理 堤坝施工 在分析开始 所有堤坝单元都不激活 在增量1分析时 施工第一层 因而在本增量开始 第一层的所有单元都被激活 然后施加第一层的自重 装配总刚矩阵 进行计算 其他各层土的填筑与此相同 通过累计每级增量分析结果可得最后计算结果 显然 结果依赖于每级填土层的厚度 layer1 layer2 layer3 1按一组增量进行分析 在某级增量或填筑开始前 确保所有拟填筑单元未激活 2在进行某级增量分析时 激活该级增量对应的填筑单元 并赋予适合于填筑过程中材料性质的本构模型 通常此时填筑材料的刚度较低 3计算由于填筑材料自重引起的节点力 并加入平衡方程右边的载荷列矢量 4组装该级增量的总刚矩阵和边界条件 解方程得该级增量位移 应变 应力 5在下级增量开始前 改变刚填筑单元的本构模型使其能表征已就位填筑材料的性质 对仅与本级填筑有关的单元节点 或新增节点 位移归零 6进行下级增量分析 数值模拟 施工过程 1随着施工的进行 单元 新填土层 逐渐被激活 并参与计算 尚未填筑的土层单元未被激活也不参与计算 总刚与载荷矢量组装时 2所有单元 填土层 都处于激活状态 只是在未被填筑时 其被赋予较小的刚度 但每一增量计算时所有单元都参与 数值模拟 网格处理的两种方式 很多软件采用第二种处理方式 layer1 layer2 layer3 Modifydisplacementsfornewfilledlayers to770 IG 0DO730IE NF NEDO730JE 1 4I ME IE JE IF IE EQ NF AND JE LT 4 GOTO720DO710J 1 IGI0 IH J IF I EQ I0 GOTO730710CONTINUE720CONTINUEIG IG 1IH IG I730CONTINUEDO770I00 1 IGI IH I00 XZ 1 X I 1 XZ 2 X I 2 CALLCYY CY0 NC0 IC0 NC1 NC NL DO740IJ 1 NC1740IC1 IJ IC NL IJ CALLCYY CY1 NC1 IC1 DZ CY1 XZ 2 DH CY1 CY0IF DH LT 0 00001 THENAA 1 0GOTO750ENDIFAA 2 0 DZ DH DZ 750CONTINUEDO760J 1 2D1 AA 1 0 DU I J 100 0U I J U I J D1UU I J UU I J D1760CONTINUE770CONTINUE 土体挖方工程 5 有限元法计算分析挖方工程 是要模拟被开挖的那部分土体去除后 地基中的应力变形状态 分析对挖坡稳定性的影响 开挖是逐步实施的 计算模拟也要分级逐步进行 一 开挖计算的主要问题 1 开挖体的处理 1 保留原来的网格体系 令挖去土体的单元的材料由土体变为空气 对空气用很小的弹性模量 2 从网格中删去已挖的土体单元 这种方法更合理 但增加了计算程序中处理的麻烦 即单元数 节点数逐级减少 ABCDEF 2 应力的释放 1 直接由挖去的各层的土重求等效的节点荷载 将其反号 成为负的荷载作用在开挖面上 但这只能看作一种粗略的估算 实际上拟开挖的面上各点在开挖过程中应力状态已经发生变化 2 对开挖工程的模拟 要计算拟开挖面上的实际应力状态 将其释放为零 而成为应力自由面 这就需要根据开挖进展不断作修正计算 当开挖到某一层时 将此时该开挖面上各点的应力分量转换为等效的节点荷载 取其负值 作用在开挖面上 以抵消开挖前该面上的应力 以此作为整个地基上的荷载 进行有限元计算 第一层 第二层 A B 由与开挖边界相邻的开挖单元的应力计算节点力 3 卸载 挖方施工时土体处于卸载状态 这里卸载是指围压降低 至于剪应力 随着开挖的进展还是增加的 仍然为加载 在选择本构模型时要注意 模型要能反映围压降低与增加对变形的不同影响 二 开挖面上应力释放的模拟 第一级开挖模拟时 必须先求得开挖前的应力 这时一般假定土体在自重应力作用下处于K0状态 即平面应变条件下为 经过第一级开挖 地面不再水平 以后各级便不能再按自重应力的方法求了 而要用有限元法计算出应力 有限元法直接求得的是单元应力 可用如下两种方法求边界应力 1由单元应力内插求节点应力 2由位移直接求节点力 在第一层开挖前 三个开挖面上的应力均为土体自重应力 三个开挖面以上土体对其下部土体所作用的等效节点力 就由自重应力求得 分别以 F I F II和 F III表示 第一阶段开挖 在开挖面上作用节点荷载 F I 去除该面以上土体 作有限元计算 得出相应的地基内应力 变形 同时也可算出第二和第三开挖面上的节点力由于第一层的开挖而产生的变化量 F II1和 F III1 F II1由紧挨在第二开挖面上方的单元劲度矩阵与第一阶段的单元节点位移增量相乘求得 此时第二开挖面上的节点力为 F II F II F II1 F III F III F III1 同理第三开挖面上的节点力为 第二阶段开挖 在第二开挖面CEFD上作用 F II 使该面成为应力自由面 从网格中删除该面以上单元 用有限元法解得位移和应力 求第二阶段开挖在第三开挖面上所产生的节点力变化 F III2 方法同前 进而求第三开挖面上当前的节点力 第三阶段开挖 就是在第三开挖面上作用节点荷载 F III F III F III F III2 F II1由紧挨在第二开挖面上方的单元劲度矩阵与第一阶段的单元节点位移增量相乘求得 M点上方紧挨的两个单元为3和4 分别用这两个单元的劲度矩阵与单元节点位移增量相乘 求得相应的节点力增量 其中只取M点的节点力增量 将这两个单元在M点的节点力增量叠加起来 就是所求的M点的 F II1 三 水下开挖 H A A H B G A H B G C F 1去掉AH上的水压力 同时施加ABGH上的水压力 2开挖第二层时 去掉B G节点由于原BG面上水压力引起的节点力 因为相关部分的土体没开挖 3BG面上B G节点以外的其他节点力可以不管 因相关单元的开挖模拟时会间接考虑 cc推求开挖释放等效结点荷载cSUBROUTINEFRUE IE NP W REAL 8RSWB 0 00DO10I 1 NP10B B XZE I 2 FLOAT NP DO40I 1 NPIP INE I J I 1IF J GT NP J 1JP 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 GOTO40DX 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 cPY 0 50 DZ SS IE 3 DX SS IE 2 W Z B PY 0 50 DZ SS IE 3 DX SS IE 2 IF PX lt 0 PX PXIF PY lt 0 PY PY A W DX DZ 12 00IF I1 NE 0 RSW I1 RSW I1 PXIF I2 NE 0 RSW I2 RSW I2 PYIF J1 NE 0 RSW J1 RSW J1 PXIF J2 NE 0 RSW J2 RSW J2 PY40CONTINUEIF NP EQ 3 GOTO60DO50IC 1 4CALLEJB IC 2 DJ DO50I 1 4L IWU INE I 2 IF DJ LT 0 DJ DJIF L GT 0 RSW L RSW L W DJ FN I IC 50CONTINUEGOTO8060CALLEBC 1 DJ DO70I 1 3L IWU INE I 2 70IF L GT 0 RSW L RSW L W DJ 3 080CONTINUERETURNENDc 空气单元 平面问题的地下结构单元 6 1 墙 混凝土墙 实体单元 钢板桩 混凝土桩或其他钢与板组合墙 梁单元 2 地锚 实体单元 弹簧单元 接触面单元 弹簧单元 弹簧单元 A B C 膜单元 膜单元 A B C 弹簧单元 A B C 弹簧单元 膜单元 膜单元 3 墙与支撑的连接 在平面应变分析中 挡土墙与支撑的连接有以下方式 Shear 0Moment 0 TiedAandBinHorizontaldirection Axialthrust 0 Solidelements TiedAandBinHorizontaldirection Beam barelements AB 简单连接 Moment 0 Axialthrust 0shear 0 TiedAandBinHorizontalandVerticaldirection Solidelements Solidelements Beam barelements AB 销连接 Axialthrust 0shear 0moment 0 Solidelements Beam barelements 完全连接 4 耦合分析中的结构单元 在耦合固结分析中 挡土墙如图a 用梁单元模拟墙如图b 由于实体单元的节点为普通节点 在固结分析中实际意味墙是透水的 水可自由流过 这与实际不符 如果要使墙为不透水 可以在墙的一边设置接触面单元如图c 如果这些接触面为非固结单元 他们将在墙的某一边提供一个不透水的隔断 abc 浅基础与桩基础的有关计算 7 面基础 浅基础 1柔性基础 在传统土力学中 分析面基础时 通常假设为两种极端情况 即理想柔性基础和理想刚性基础 在数值计算时如果也用这种假设 则在计算过程中可以不包括基础本身 而以边界条件形式考虑 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 8m 塔顶离中心线已达5 27m 倾斜5 5 成为危险建筑 1990年1月4日被封闭 除加固塔身外 用压重法和取土法进行地基处理 为了抢救斜塔 1930年曾给塔基灌浆 做加固处理 见效甚微 1965年和1973年意大利政府曾出高价1万至2 5万英镑向世界各国征集纠斜方案 1990年意大利政府为了保证游人的安全 关闭了斜塔 此后继续征集方案 最终还是由聪明的中国人的方案中标施工 采取了地基挖出100吨土等措施处理了地基 经测定塔再未倾斜 还扶正了5毫米 因此近年来斜塔又向游人开放 可登塔观光 掏土过程模拟 将单元6的刚度减小至0 在单元6的上 下边逐渐施加等值反向的垂直节点力 直到单元体积减小约5 然后回复单元刚度 连续对单元7 8 9 10和11重复上述过程 以模拟陶土器的插入过程 每一步都会使塔的倾斜减小 当开挖单元12时 塔的倾斜反而增加 即临界线右边的单元开挖导致负响应 因此 开挖单元11后 分析应重新开始 通过连续开挖单元10 9 8 7和6 模拟掏土器的撤回 每一步 塔的响应为正 倾斜减小 重复陶土器的插入和撤回过程 再次确认单元12开挖导致负响应 倾斜变化 deg 基础沉降 mm 6789101110987678910111098767891011 300 200 100 0 1 0 2 0 3 0 4 Retraction Insetion south north north 5桩基础 单桩 在桩附近用较小尺度的单元 或使用接触面单元 桩体用实体单元 仅受垂直荷载作用 可不用接触面单元 受侧向荷载或组合荷载作用 有必要用接触面单元 垂直荷载 Nointerfaceelement Ks Kn 105kN m3 Ks Kn 103kN m3 8000600040002000 00 020 040 060 080 10 120 14 轴对称或平面对称 水平荷载 Nointerfaceelement Withinterfaceelement 4000300020001000 00 020 040 060 080 1 kN m 三维分析 05101520 05101520 Nogapping gapping Notension tension 800 4000400 水平应力kPa 桩深m 桩深m 0 0200 020 040 060 08 水平位移m 群桩 叠加技术 刚性桩帽 在原理上 群桩基础可用三维有限元分析 但计算成本过大 因此 必须对实际问题予以简化