基坑开挖数值模拟.docx

1、基坑开挖数值模拟7 数值模拟7.1数值模拟方法简介数值模拟技术作为一种研究手段， 已经被广泛的应用于各行各业 领域的研究中。目前，数值分析方法主要分为二大类：一类是以有限 差分法为代表，其特点是直接求解基本方程和相应的定解条件的近似 解；另一类数值分析方法是首先建立和原问题基本方程及相应定解条 件等效的积分方法，然后据之建立近似解法。LS-D YNA乍为世界上最著名的通用显示动力分析程序，能够模拟 真实世界的各种复杂问题， 特别适合求解各种二维三维非线性结构的 高速碰撞， 爆炸和金属成型等非线性动力冲击问题， 同时可以求解传 热，流体及流固耦合问题， 在工程应用如汽车安全设计，武器系统设 计，

2、金属成型，跌落仿真等领域被广泛应用。 本次采用 ANSYS/LS-DYN，A 进行混凝土支撑梁结构爆破拆除数值模拟研究。在ANSYS/LS-DYN环 境下，数值模拟的实现总体上分为两个过程： 在ANS丫芽建立结构实 体模型，完成有限元网格的划分， 输出有限元模型信息即输出关键字 文件；编辑关键字文件，在 DYNA环境下完成对结构倒塌过程的数值 模拟计算。对结构有限元模型的建立过程， 数值模拟中采用的钢筋和混凝土 材料模型、接触方式等各种计算控制项进行了阐述。LS-D YNA程序中主要提供如下几种计算方法：（ 1） Lagrange 算法坐标固定在物质上或者说随物质一起运动和变形， 处理自由面和

3、 物质界面非常直观， 由于网格始终对应物质， 因此能够精确的跟踪材 料边界和描述物质之间的界面，这是 Lagrange 算法的主要优点。但 是，由于网格随材料流动而变形，一旦网格变形严重，就会引起数值 计算的不稳定， 甚至使得计算无法继续进行 （如发生负体积或复杂声 速等问题）。因此， Lagrange 算法在处理大变形大位移问题时，有其 无法克服的弊端。（ 2）Euler 算法 网格被固定在空间，是不变形的。物质通过网格边界流进流出， 物质的大变形不直接影响时间步长的计算。 因此，欧拉算法在处理大 变形问题方面具有优势。欧拉方法通过输运项计算体积、质量、动量 和能量的流动。 欧拉计算可以直接

4、通过在离散化格式中包括迁移导数 项进行， 或通过二步操作完成。 二步法操作的第一步主要是拉格朗日 计算，第二步输运阶段是重分计算网格相当于回到它的原来状态。LS-D YNA程序采用后一种方法。欧拉算法的缺点是网格中物质边界不 清晰，难以捕捉各物质界面。（3）ALE方法吸取了欧拉法和拉格朗日法两种方法的优点。 ALE算法能够进行 自动重分网格操作。它包括拉格朗日时间步，然后是一个输运步。输 运步可以采用三种方法： 1. 发生合理的网格变形时空间网格不再重 分（拉格朗日）；2. 发生严重的网格变形时重分成原始形状 （欧拉）； 3. 发生严重的网格变形时重分为合理的形状， 因此允许网格拓扑（拉 格朗

5、日和欧拉）。混凝土是土木工程结构中应用极为广泛的材料， 其最本质的特点 是材料组成的不均匀性，并且存在初始微裂缝。从混凝土受单轴压力 时的应力应变关系来看，混凝土卸载时有残余变形，不符合弹性关系； 如果对其应用弹塑性本构关系，又很难精确定义屈服条件。此外，混 凝土在到达应力顶峰后，其应力-应变关系曲线有一下降段，即存在 应变软化现象，所有这些都给建立混凝土的本构关系带来困难。 多年以来，众多学者进行了大量的试验和理论研究， 提出了各种各样的混 凝土本构模型。第一：混凝土本构模型可以分为下面几种：（1） 线弹性类本构模型。线弹性类本构模型是以弹性力学为基 础的模型，当混凝土无裂缝时，将混凝土看成

6、线弹性匀质材料而采用 线弹性本构模型。虽然混凝土的变形特性是非线性的， 但在一些特定 的情况下（比如描述混凝土受拉时的工作性能），采用线弹性类本构模 型进行分析还是有足够的精度的，其线弹性本构关系可用广义虎克定 律来表示：ij - Cijkl ；kl （ 1）式中，Cijk 1为材料弹性常数，为四阶张量，共有 81个常数。按照 材料假设的不同，又可分为各向异性本构模型、正交各向异性本构模 型、各向同性本构模型等，其中Cjki根据材料的不同而变换。（2） 塑性理论类本构模型：塑性理论类本构模型是以塑性流动 理论为基础，代表性的模型主要有： Mises条件的模型、理想弹塑性脆性断裂模型、 应变或工

7、作硬化塑性理论模型等。 模型中考虑了混凝 土加载路径和混凝土的硬化，在混凝土的应力 - 应变全曲线中，有上 升段和下降段。自从Drucke公设和nyushin公设出现之后，经典塑 性力学得到飞速发展， 混凝土塑性力学模型也是基于这些公设建立的。 以塑性理论为基础的混凝土本构模型， 在对其加载面， 包括初始屈服 面，后续加载面和破坏包络面等特征面的研究中， 这些特征面若以应 力空间来表示时，当应力达到屈服后，材料发生应力松弛；若以应变 空间表示时，当应变达到松弛面后，材料发生应变松弛。基于应力状 态屈服面或破坏包络面的塑性理论类型的本构模型有弹性 - 全塑性模 型、线弹性 -硬化塑性 -断裂模型

8、等； 基于松弛面的塑性理论类型的本 构模型有塑性模型、塑性断裂模型、硬化断裂模型等。所有这些模型 所做的假设与混凝土的实际性能还存在很大的差别， 而且模型的表达 式和计算均较复杂，目前还不便于应用。( 3)其他力学理论类本构模型，许多学者还以新型交叉的力学 分支的理论为基础研究混凝土的本构模型。内时理论模型：内时理论模型最初由 Valanis 于 1971 年提出， 其基本概念为：塑性和粘塑性材料内任一点的现时应力状态是该点整 个变形领域内和温度历史的泛函， 而特别重要的是该历史是用一个取 决于变形中的材料特性和变形程度的内时 (Intrinsic Time) 来量度 的。这种模型采用了非弹性

9、应变能逐渐积累的方法而不需考虑塑性理 论中的屈服面和流动法则， 所以该理论尤其适合没有屈服面的混凝土 材料。由于内时理论能描述混凝土的复杂变形的历史， 因而为各国学 者所重视。但由于表达式过多，确定参数又不容易，所以对其推广和 应用仍有大量工作要做。损伤理论模型：针对混凝土材料本身是一种具有固有缺陷微裂 纹的特点， 很多学者将损伤力学引用到混凝土本构模型的建立中。 损 伤力学研究材料或构件从原生缺陷到形成客观微裂纹直至断裂的过 程。也就是通常指的微裂纹的萌生、 扩展或演变、宏观微裂纹的形成、 裂纹的稳定扩展和失稳扩展全过程。 损伤力学主要是在连续介质力学 和热力学的基础上， 用固体力学方法研究

10、材料的宏观力学性能的演变 直到破坏的全过程。 20世纪 70年代末期，损伤力学局限于研究材料 在宏观裂纹出现以前的阶段， 当宏观裂纹出现以后， 则采用断裂力学 的理论和方法进行研究，这是无耦合的分析方法。实际上，在宏观裂 纹出现后，材料的损伤对裂纹尖端附近及其它区域的应力和应变均有 影响。因此，合理的方法应该将损伤耦合到本构方程中进行分析和计 算。这样由于本构方程中将有关的力学参数和损伤进行了耦合， 所以 分析和计算就变得更为复杂。非线弹性类本构模型： 为了克服线弹性模型的缺点， 然后产生了 以割线模型形式的非线性弹性类本构模型。 这类模型中， 具有代表性 的是超弹性模型和亚弹性模型。除了上述

11、模型外， 还有以许多种理论相结合到处的混凝土的本构 模型。第二：LS-D YNA程序中提供的混凝土材料模型：在材料模型方面，LS-DYNA目前拥有近150余种金属和非金属材料模型，涵盖了弹性、弹塑性、超弹、泡沫、玻璃、土壤、混凝土、流体、炸药等材料模型以及多种状态方程，可以考虑材料的失效、损 伤、粘性、蠕变、与温度相关、与应变率相关等材料性质。此外，程 序还支持用户自定义材料功能。LS-D YNA程序中提供的混凝土材料模型主要有以下几种：塑性随动模型 *MAT-PLASTIC-KINEMATICH_J_C模型、*MAT_DRUCKER_PRAGER 模型和 *MAT_SOIL_AND_FO模型

12、。本论文有限元模型中的混凝土本构模型选用 LS-D YNA中的塑性随动模型 *MAT-PLASTIC-KINEMATQ塑性随动模型（*MAT_PLASTIC_KINEMATIC这是一种各向同性、随动硬化或各向同性和随动硬化混合模型， 与应变率相关， 可考虑实 效，参数简单，较容易确定。7.2钢筋混凝土支撑梁结构爆破拆除有限元模型的建立7.2.1 爆破拆除数值模拟基本假设考虑到建立整体钢筋混凝土支撑梁结构模型过程的复杂性， 计算 结果的可行性及准确性， 数值模拟中对模型进行了简化， 模型的简化 出于以下基本假设：（ 1）不考虑混凝土支撑梁内的炸药爆炸对整体结构的作用。如 果模拟中按照实际情况建立

13、模型， 单元数目巨大， 受到当前计算机能 力的限制，无法实现数值计算。由于爆炸对钢筋梁的变形影响较小， 所以在模拟中可以直接删除爆破药包附近的混凝土， 并没有考虑钢筋梁结构在炸药爆炸下的响应问题（2）钢筋混凝土梁在爆破后，剩余部分的表面是光滑的，由于 上面假设是直接删除爆破部分，所以无法模拟混凝土爆破后的形状。（3）钢筋混凝土的支柱地面没有相对位移，简化为完全固结， 即不考虑钢筋梁的坍塌与支撑柱和地面之间的相互作用。7.2.2钢筋混凝土有限元模型建立过程采用ANSYS建立有限元模型时，采用的数值单位均为国际单位， 即长度单位为m时间单位为s,质量的为kg,压力的单位为Pa,速 度单位为 m/

14、s 。为了由实际情况对比实验模拟的准确情况, 我们选用工程实例作 为有限元模拟的对象。 综合考虑有限元计算的精确性和时效性, 来确 定单元尺寸和单元数量。7.2.3 钢筋混凝土模型的选取采用数值模拟方法进行钢筋混凝土结构爆破拆除分析前, 首先要 就研究的具体问题选择好用于模拟钢筋混凝土材料的物理模型。 当前, 数值计算中, 用于模拟混凝土钢筋材料的有限元模型主要有两类： 第 一种是把钢筋和混凝土单元各自划分为足够小的单元, 分别考虑钢筋 和混凝土的贡献, 称为分离式模型。 第二种也是把钢筋和混凝土单元 包含在一个单元之中, 统一考虑钢筋和混凝土的作用, 称为整体式模（ 1）整体式模型在整体式模

15、型中， 将钢筋弥散于整个单元中， 并把单元视为连续 均匀的材料。 钢筋对整个结构的贡献， 可以通过调整单元的材料力学 性能参数来体现，例如提高材料的屈服强度、材料的弹性模量等。其 优点是建模方便， 分析效率高， 但是缺点是不适用于钢筋分布较不均 匀的区域， 且得到钢筋内力状态比较困难。 主要用于有大量钢筋且钢 筋分布较均匀的构件中。（2） 分离式模型，位移协调利用空间梁单元beam161建立钢筋模型，和混凝土单元共用节点。 其优点是建模方便， 可以任意布置钢筋并可直观获得钢筋的内力。 缺 点是建模比整体式模型要复杂， 需要考虑共用节点的位置， 且容易出 现应力集中问题。（ 3）分离式模型，界面

16、单元 前两种混凝土和钢筋组合方法假设钢筋和混凝土之间位移完全 协调，没有考虑钢筋和混凝土之间的位移， 而通过加入界面单元的方 法，可以进一步提高分析的精度。同样利用空间梁单元beam161建立 钢筋模型。不同的是混凝土单元和钢筋单元之间利用弹簧模型来建立 连接。不过一般钢筋混凝土结构中钢筋和混凝土之间都有比较良好的 锚固，钢筋和混凝土之间滑移带来的问题不是很严重， 一般不必考虑。本文，在爆破拆除建筑结构的数值模拟中， 兼顾模型的真实性和 建模的方便性， 数值模拟中采用分离式钢筋混凝土模型， 钢筋和混凝 土之间共节点， 不考虑它们之间的滑动。 并对结构中钢筋的分布位置及数量进行了简化 7.2.4

17、 钢筋混凝土的单元类型LS-D YNA具有丰富的单元库，具有二维、三维实体单元，薄、厚 壳单元以及 A LE、Eulirian 和 Lagrangian 单元等，各类单元又有多 种算法可供选择，具有大位移、大应变和大转动功能，单元积分采用 沙漏粘性阻尼以克服零能模式。模拟中，混凝土单元选择 SOLID164 单元，钢筋单元选择 BEAM16单元。SOLID164单元，为8节点三维实体单元，每个节点具有 9个自 由度，即 UX UY uz vx vy vz AX AY AZ但是只有位移是 实际意义上的自由度。 该单元提供了单点积分算法和全积分单元算法 其默认为单点积分算法。 该单元可采用 Lag

18、range 列式 也可采用 ALE 列式 当采用 ALE 列式时 在关键字中需要定义单元算法关键字 *SECTION_SOLID_ALE单元支持大部分的 LS-D YNA材料算法。BEAM16单元，为二维梁单元，用三个节点定义，I和J确定梁 的轴向， K 确定横截面主轴方位。具有 Hughes-Liu 及 B-S 算法，其 默认算法为 Hughes-Liu 算法。单元可以定义不同的积分算法，默认 值为2x2Guass积分。对于Hughes-Liu算法选项，如果截面积分规 则 Int.rule-arbitry.section 取默认值 0，则程序会采用梁单元跨 中的一组积分点来模拟矩形及圆形的截

19、面。725材料失效方式的控制在有限元数值模拟中，由于程序本身是基于连续介质力学理论基 础上形成的，因此，材料的失效是靠删除单元来实现的。单元达到屈 服强度后，就被从模型中删除，不再参与计算，所以，对于材料失效 强度的控制决定着数值模拟计算的准确性。材料失效现象比较复杂，因此，出现了各种失效理论模型，目前, 比较成熟的强度理论模型主要有以下四种：第一强度理论模型，即最大拉应力理论模型;G 一 丿 （3-1）第二强度理论模型，即最大伸长线应变理论模型;(3-2)第三强度理论模型，即最大剪应力理论模型;第四强度理论模型，即形状改变比能理论模型（Von Mise应力）(3-4)四式中，J L许用应力，

20、屈服极限二s/安全系数；J，二2，二3-三个主应力。对于铸铁、石料、混凝土和玻璃等脆性材料，通常以断裂形式失 效，可采用第一和第二强度理论；对于碳钢、铜、铝等塑性材料，通 常以屈服形式失效，可以采用第三和第四强度理论。根据以上论叙，对于混凝土材料，论文中采用材料控制项中 FS 和额外控制项*MAT_ADD_EROSION合控制材料的失效。对于钢筋材 料则采用材料模型*MAT_PLASTIC_KINEMAT中提供的FS参数控制。7.2.6 接触方式的选择LS-D YNA提供丰富的接触方式可供选择，可以求解各种柔性体之 间、刚性体之间以及柔性体和刚性体之间的接触问题， 并可分析接触 表面的静动摩擦

21、力、固连失效以及流体与固体的界面等问题。程序中处理动态的接触 -碰撞界面主要采用如下三种算法。（1）动力约束法 动力约束法是最早采用的接触算法， 该方法在每一步修正结构的 构形之前，判断每一个没有与主表面接触的从节点是否会在该时步内 贯穿主表面。如有从节点贯穿主表面，则将时步缩小，使那些从节点 都不贯穿主表面，而其中有的刚到达主表面，在下一时步开始，对刚 到达主表面的从节点施加动力约束条件。 对所有已经与主表面接触的 从节点都施加约束条件，以保持从节点与主表面的接触。此外，检查 与主表面接触的从节点所属的单元是否存在受拉交界面力。 如果受拉 面力，则用释放条件使从节点脱离主表面。由于此算法比较

22、复杂，后 来只用于固连接触问题，即没有释放条件。（2）分配参数法该方法仅用于滑动界面接触问题的处理， 其原理是：将所有正在 接触的从单元质量的一半分配到主表面上， 同时由每个单元的应力确 定作用在接受质量的主表面面积上的分布压力。 在完成质量和压力的 分配后，程序修正主表面的加速度。 然后对从节点的加速度和速度施 加约束， 以保证从节点沿主表面运动。 这一算法不允许从节点穿透主 表面，适合相对滑动但没有分离的表面之间的接触分析。（ 3）对称罚函数法对称罚函数法为LS-DYNA程序的缺省接触算法，该方法的基本原 理为：每一时步先检查各从节点是否穿透主表面， 没有穿透则对该节 点不作任何处理。 如

23、果穿透， 则在该从节点与被穿透主表面之间引入 一个较大的界面接触力，其大小与穿透深度、主片刚度成正比，称为 罚函数值。它的物理意义相当于在从节点和被穿透的主表面之间放置 一个法向弹簧， 以限制从节点对主表面的穿透。 所谓对称罚函数法是 指程序在对所有的主节点按如上步骤处理， 其算法与从节点一样。 对 称罚函数法，很少激起网格的零能模式，没有噪声，这是由于算法具 有对称性、动量守恒准确， 不需要碰撞和释放条件。罚函数值大小受 到稳定性限制。 若计算中发现明显穿透， 可放大罚函数值或缩小时间 步长来调节。在LS-DYNA程序中，接触类型可以分为三类，即单面接触、节点- 表面接触和表面 - 表面接触

24、。（1）单面接触单面接触包括的接触类型有 SS ASSC AG ASS2D ESS单面 接触可以用于一个物体表面各部分的自相接触或与另一个物体的表 面接触，在单面接触中，程序会自动判定模型中发生表面接触的位置。 因此，单元接触的定义是最简单的，无需定义源面和目标面，当定义 好单面接触时， 它允许一个模型的所有外表面都可能接触， 有利于处 理预先不知道接触表面的自身接触或大变形问题。（2）点- 面接触LS-D YNA程序中的点-面接触类型有 NTS ANTS RNTRTDNSTNTSENTSDRAWBEADNTS等。对于点-面接触，接触节点将穿透目标面， 这种接触类型用于一般情况下的两个物体表面

25、之间的接触。点 -面接 触需要定义 target 面和 contact 面，原则如下：平面或凹面为 target 面，凸面为 contact 面；网格较粗的面作为 target 面，网格较细为 contact 面。（3）面-面接触LS-DYNAg序中的面-面接触类型包括 STS OSTS ASTS ROTR TDSS TSTS ESTS SE FSTS FOSS TSES 等。当一个物体的表面 穿透另一个物体表面时需用面 -面接触，这种接触是最常用的接触类 型，常用于任意形状且存在较大接触面积的物体之间的接触问题。 该 接触类型对于物体间有大量相对滑移时很有效。接触碰撞问题属于最困难的非线性问

26、题之一， 因为在接触碰 撞问题中的响应是不平滑的。 当发生碰撞时， 垂直于接触界面的速度 是瞬时不连续的。 由于分析本身非常复杂， 接触情况变化多端。 因此， 计算中选用了 LS-DYNA提供的CONTACT-ERODING-SINGLE-F接触计 算模型。该接触计算模型可以自动搜索接触面，判断接触，并可以处理侵蚀、断裂等复杂边界变化情况。接触刚度选用 1.1 ，材料的摩擦 系数统一设定为 0.6 。7.2.7 重要参数的设定对于梁结构的每个端面和支撑柱的地面部分以及模拟地面的地 板添加全方位的面约束。在生成关键字后，在关键字中改为 Z 方向的，重力加速度为 9.8m/s 。计算0到5秒内的爆

27、破过程，将 DATABASE_BINARY_D3PI设置 为 0.025s,将 DATABASE_BINARY_D3TI也设置为 0.025s。由于计算 机的计算能力限制每秒钟不能超过 1000个 plot 文件，也就是说时间 不能低于 0.005s。728钢筋混凝土支撑梁结构爆破拆除数值模拟研究流程图图7.1钢筋混凝土支撑梁结构爆破拆除数值模拟研究流程图图7.1为本次数值模拟研究的流程图，整个建模以及计算过程都是在ANSYS12.0F完成的，建立的数值计算模型如图 7.2所示。7.3.有限元模拟计算结果以及分析7.3.1爆破过程分析为实现对爆破设计中的延时起爆顺序的模拟， 将支撑梁的不同爆区

28、建立不同的部件，每个部件通过*MAT_ADD_EROSION腔制混凝土失效的时间从而来模拟 爆破过程。每个爆区的延时时间为 25ms,整个延时起爆过程如图 7.3图7.20LB-DYIMine* fenpMt所示，La-DYhKufflWkiMJl TtM LB-DYM mi- QMtnaLB-DYIM !HH bnp/ m- a .1LB-DYMiHFFlnpJ: nrw Qnn“OYMAszinpwi rw omL8*0YMAuie Input rw- o；Lt-OYHAme* Input rw- 02ML8*OYMAuie0YNA cser Inputtw- DSM90LS-OYNAme

29、r Input tw- e*L9-0YM4 uwc mput rm* sr&LS-OYMA seWWLS-OYMAww InputLS-OYHAmer Inputrw- 2C2732混凝土单元的最大主应力随时间的变化此单元位于第一部分内，所以一开始就有一个压力突变，此单元为掉落在地上，是与钢筋连在一起，所以一直会有压力的变化。在0.4s时达到最大值。Element no.A. 9129(g+山二恥S切 匕仞ffidoEudErlExEW此单元为第二部分的，在一开始会有强烈的压力突变，因为在中Time(s)图表4-1间部分，第二峰值是因为此单元掉落在地上， 之后的波动则是单元之间的碰撞。第三部分

30、，此单元紧挨立柱，所以一开始并无突变，在 0.7s左右达到最大值，是因为单元一部分位于钢筋分开以至于掉落碰撞到住上，之后平稳。Time(s)图表4-18此单元在第四部分内，而属于中间部分,=LS-DYNA user input ee0.80.60.40.20gTime(s)图表4-19的压力突变在0.7s左右有强烈的突变是因为此单元掉落地上，之后 的突变是因为她压在一个单元上面。A 18143(g+山二nlssul巴scbbh 匚一UJnE 一第五部分内，此单元最后落在在支柱上面，且距离支柱的顶部有 很短的距离，在一开始之后又两个突变的峰值，第一个是删除后的变 形，第二是碰撞到支柱顶部的，之后

31、不稳定是因为此单元还是有一部 分连接着钢筋，切支柱是多个梁的交点。第六部分内，此单元与前边的分布距离较近，所以有两个相邻的峰值一个是前面部分删除引起的，一个是本身引起的。在后在 0.8s左右的突变是掉落。启 22766山)fficltt孙 mhsrtd.o匚 UCLEnEx如左图表4-21LS-DYNA user input ele-.J11 * 一I 1， T-W 一L L AF E ElefnentnoA 26653Timet s)(守山二Ed茶常怎sduu匸亠EnEM毎壬第七部分内，峰值是因为变形而产生的，没有一开始发生突变是因为此部分位于边沿，没有直接接触，且距离支柱较近，此面不是碰 击地面的那部分，并且有一部分连接到钢筋上。第八部