基坑开挖数值模拟.docx

7数值模拟7.1数值模拟方法简介数值模拟技术作为一种研究手段，已经被广泛的应用于各行各业领域的研究中。目前，数值分析方法主要分为二大类：一类是以有限差分法为代表，其特点是直接求解基本方程和相应的定解条件的近似解；另一类数值分析方法是首先建立和原问题基本方程及相应定解条件等效的积分方法，然后据之建立近似解法。LS-DYNA作为世界上最著名的通用显示动力分析程序，能够模拟真实世界的各种复杂问题，特别适合求解各种二维三维非线性结构的高速碰撞，爆炸和金属成型等非线性动力冲击问题，同时可以求解传热，流体及流固耦合问题，在工程应用如汽车安全设计，武器系统设计，金属成型，跌落仿真等领域被广泛应用。本次采用ANSYS/LS-DYNA，进行混凝土支撑梁结构爆破拆除数值模拟研究。在ANSYS/LS-DYNA环境下，数值模拟的实现总体上分为两个过程：在ANSYS中建立结构实体模型，完成有限元网格的划分，输出有限元模型信息即输出关键字文件；编辑关键字文件，在DYNA环境下完成对结构倒塌过程的数值模拟计算。对结构有限元模型的建立过程，数值模拟中采用的钢筋和混凝土材料模型、接触方式等各种计算控制项进行了阐述。LS-DYNA程序中主要提供如下几种计算方法：（1）Lagrange算法坐标固定在物质上或者说随物质一起运动和变形，处理自由面和物质界面非常直观，由于网格始终对应物质，因此能够精确的跟踪材料边界和描述物质之间的界面，这是Lagrange算法的主要优点。但是，由于网格随材料流动而变形，一旦网格变形严重，就会引起数值计算的不稳定，甚至使得计算无法继续进行（如发生负体积或复杂声速等问题）。因此，Lagrange算法在处理大变形大位移问题时，有其无法克服的弊端。（2）Euler算法网格被固定在空间，是不变形的。物质通过网格边界流进流出，物质的大变形不直接影响时间步长的计算。因此，欧拉算法在处理大变形问题方面具有优势。欧拉方法通过输运项计算体积、质量、动量和能量的流动。欧拉计算可以直接通过在离散化格式中包括迁移导数项进行，或通过二步操作完成。二步法操作的第一步主要是拉格朗日计算，第二步输运阶段是重分计算网格相当于回到它的原来状态。LS-DYNA程序采用后一种方法。欧拉算法的缺点是网格中物质边界不清晰，难以捕捉各物质界面。（3）ALE方法吸取了欧拉法和拉格朗日法两种方法的优点。ALE算法能够进行自动重分网格操作。它包括拉格朗日时间步，然后是一个输运步。输运步可以采用三种方法：1. 发生合理的网格变形时空间网格不再重分（拉格朗日）；2. 发生严重的网格变形时重分成原始形状（欧拉）；3. 发生严重的网格变形时重分为合理的形状，因此允许网格拓扑（拉格朗日和欧拉）。混凝土是土木工程结构中应用极为广泛的材料，其最本质的特点是材料组成的不均匀性，并且存在初始微裂缝。从混凝土受单轴压力时的应力应变关系来看，混凝土卸载时有残余变形，不符合弹性关系；如果对其应用弹塑性本构关系，又很难精确定义屈服条件。此外，混凝土在到达应力顶峰后，其应力-应变关系曲线有一下降段，即存在应变软化现象，所有这些都给建立混凝土的本构关系带来困难。多年以来，众多学者进行了大量的试验和理论研究，提出了各种各样的混凝土本构模型。第一：混凝土本构模型可以分为下面几种：（1）线弹性类本构模型。线弹性类本构模型是以弹性力学为基础的模型，当混凝土无裂缝时，将混凝土看成线弹性匀质材料而采用线弹性本构模型。虽然混凝土的变形特性是非线性的，但在一些特定的情况下(比如描述混凝土受拉时的工作性能)，采用线弹性类本构模型进行分析还是有足够的精度的，其线弹性本构关系可用广义虎克定律来表示：（1）式中，为材料弹性常数，为四阶张量，共有81个常数。按照材料假设的不同，又可分为各向异性本构模型、正交各向异性本构模型、各向同性本构模型等，其中根据材料的不同而变换。（2）塑性理论类本构模型：塑性理论类本构模型是以塑性流动理论为基础，代表性的模型主要有：Mises条件的模型、理想弹塑性脆性断裂模型、应变或工作硬化塑性理论模型等。模型中考虑了混凝土加载路径和混凝土的硬化，在混凝土的应力-应变全曲线中，有上升段和下降段。自从Drucke公设和yushin公设出现之后，经典塑性力学得到飞速发展，混凝土塑性力学模型也是基于这些公设建立的。以塑性理论为基础的混凝土本构模型，在对其加载面，包括初始屈服面，后续加载面和破坏包络面等特征面的研究中，这些特征面若以应力空间来表示时，当应力达到屈服后，材料发生应力松弛；若以应变空间表示时，当应变达到松弛面后，材料发生应变松弛。基于应力状态屈服面或破坏包络面的塑性理论类型的本构模型有弹性-全塑性模型、线弹性-硬化塑性-断裂模型等；基于松弛面的塑性理论类型的本构模型有塑性模型、塑性断裂模型、硬化断裂模型等。所有这些模型所做的假设与混凝土的实际性能还存在很大的差别，而且模型的表达式和计算均较复杂，目前还不便于应用。（3）其他力学理论类本构模型，许多学者还以新型交叉的力学分支的理论为基础研究混凝土的本构模型。内时理论模型：内时理论模型最初由Valanis于1971年提出，其基本概念为：塑性和粘塑性材料内任一点的现时应力状态是该点整个变形领域内和温度历史的泛函，而特别重要的是该历史是用一个取决于变形中的材料特性和变形程度的内时(Intrinsic Time)来量度的。这种模型采用了非弹性应变能逐渐积累的方法而不需考虑塑性理论中的屈服面和流动法则，所以该理论尤其适合没有屈服面的混凝土材料。由于内时理论能描述混凝土的复杂变形的历史，因而为各国学者所重视。但由于表达式过多，确定参数又不容易，所以对其推广和应用仍有大量工作要做。损伤理论模型：针对混凝土材料本身是一种具有固有缺陷微裂纹的特点，很多学者将损伤力学引用到混凝土本构模型的建立中。损伤力学研究材料或构件从原生缺陷到形成客观微裂纹直至断裂的过程。也就是通常指的微裂纹的萌生、扩展或演变、宏观微裂纹的形成、裂纹的稳定扩展和失稳扩展全过程。损伤力学主要是在连续介质力学和热力学的基础上，用固体力学方法研究材料的宏观力学性能的演变直到破坏的全过程。20世纪70年代末期，损伤力学局限于研究材料在宏观裂纹出现以前的阶段，当宏观裂纹出现以后，则采用断裂力学的理论和方法进行研究，这是无耦合的分析方法。实际上，在宏观裂纹出现后，材料的损伤对裂纹尖端附近及其它区域的应力和应变均有影响。因此，合理的方法应该将损伤耦合到本构方程中进行分析和计算。这样由于本构方程中将有关的力学参数和损伤进行了耦合，所以分析和计算就变得更为复杂。非线弹性类本构模型：为了克服线弹性模型的缺点，然后产生了以割线模型形式的非线性弹性类本构模型。这类模型中，具有代表性的是超弹性模型和亚弹性模型。除了上述模型外，还有以许多种理论相结合到处的混凝土的本构模型。第二：LS-DYNA程序中提供的混凝土材料模型：在材料模型方面，LS-DYNA目前拥有近150余种金属和非金属材料模型，涵盖了弹性、弹塑性、超弹、泡沫、玻璃、土壤、混凝土、流体、炸药等材料模型以及多种状态方程，可以考虑材料的失效、损伤、粘性、蠕变、与温度相关、与应变率相关等材料性质。此外，程序还支持用户自定义材料功能。LS-DYNA程序中提供的混凝土材料模型主要有以下几种：塑性随动模型*MAT-PLASTIC-KINEMATIC、H_J_C模型、*MAT_DRUCKER_PRAGER模型和*MAT_SOIL_AND_FOAM模型。本论文有限元模型中的混凝土本构模型选用LS-DYNA中的塑性随动模型*MAT-PLASTIC-KINEMATIC。塑性随动模型（*MAT_PLASTIC_KINEMATIC）。这是一种各向同性、随动硬化或各向同性和随动硬化混合模型，与应变率相关，可考虑实效，参数简单，较容易确定。7.2 钢筋混凝土支撑梁结构爆破拆除有限元模型的建立7.2.1爆破拆除数值模拟基本假设考虑到建立整体钢筋混凝土支撑梁结构模型过程的复杂性，计算结果的可行性及准确性，数值模拟中对模型进行了简化，模型的简化出于以下基本假设：（1）不考虑混凝土支撑梁内的炸药爆炸对整体结构的作用。如果模拟中按照实际情况建立模型，单元数目巨大，受到当前计算机能力的限制，无法实现数值计算。由于爆炸对钢筋梁的变形影响较小，所以在模拟中可以直接删除爆破药包附近的混凝土，并没有考虑钢筋梁结构在炸药爆炸下的响应问题。（2）钢筋混凝土梁在爆破后，剩余部分的表面是光滑的，由于上面假设是直接删除爆破部分，所以无法模拟混凝土爆破后的形状。（3）钢筋混凝土的支柱地面没有相对位移，简化为完全固结，即不考虑钢筋梁的坍塌与支撑柱和地面之间的相互作用。7.2.2钢筋混凝土有限元模型建立过程采用ANSYS建立有限元模型时，采用的数值单位均为国际单位，即长度单位为m，时间单位为s，质量的为kg，压力的单位为Pa，速度单位为m/ s。为了由实际情况对比实验模拟的准确情况，我们选用工程实例作为有限元模拟的对象。综合考虑有限元计算的精确性和时效性，来确定单元尺寸和单元数量。7.2.3钢筋混凝土模型的选取采用数值模拟方法进行钢筋混凝土结构爆破拆除分析前，首先要就研究的具体问题选择好用于模拟钢筋混凝土材料的物理模型。当前，数值计算中，用于模拟混凝土钢筋材料的有限元模型主要有两类：第一种是把钢筋和混凝土单元各自划分为足够小的单元，分别考虑钢筋和混凝土的贡献，称为分离式模型。第二种也是把钢筋和混凝土单元包含在一个单元之中，统一考虑钢筋和混凝土的作用，称为整体式模型。（1）整体式模型在整体式模型中，将钢筋弥散于整个单元中，并把单元视为连续均匀的材料。钢筋对整个结构的贡献，可以通过调整单元的材料力学性能参数来体现，例如提高材料的屈服强度、材料的弹性模量等。其优点是建模方便，分析效率高，但是缺点是不适用于钢筋分布较不均匀的区域，且得到钢筋内力状态比较困难。主要用于有大量钢筋且钢筋分布较均匀的构件中。(2)分离式模型，位移协调利用空间梁单元beam161建立钢筋模型，和混凝土单元共用节点。其优点是建模方便，可以任意布置钢筋并可直观获得钢筋的内力。缺点是建模比整体式模型要复杂，需要考虑共用节点的位置，且容易出现应力集中问题。（3）分离式模型，界面单元前两种混凝土和钢筋组合方法假设钢筋和混凝土之间位移完全协调，没有考虑钢筋和混凝土之间的位移，而通过加入界面单元的方法，可以进一步提高分析的精度。同样利用空间梁单元beam161建立钢筋模型。不同的是混凝土单元和钢筋单元之间利用弹簧模型来建立连接。不过一般钢筋混凝土结构中钢筋和混凝土之间都有比较良好的锚固，钢筋和混凝土之间滑移带来的问题不是很严重，一般不必考虑。本文，在爆破拆除建筑结构的数值模拟中，兼顾模型的真实性和建模的方便性，数值模拟中采用分离式钢筋混凝土模型，钢筋和混凝土之间共节点，不考虑它们之间的滑动。并对结构中钢筋的分布位置及数量进行了简化。7.2.4钢筋混凝土的单元类型LS-DYNA具有丰富的单元库，具有二维、三维实体单元，薄、厚壳单元以及ALE、Eulirian和Lagrangian单元等，各类单元又有多种算法可供选择，具有大位移、大应变和大转动功能，单元积分采用沙漏粘性阻尼以克服零能模式。模拟中，混凝土单元选择SOLID164单元，钢筋单元选择BEAM161单元。SOLID164单元，为8节点三维实体单元，每个节点具有9个自由度，即UX，UY，UZ，VX，VY，VZ，AX，AY，AZ，但是只有位移是实际意义上的自由度。该单元提供了单点积分算法和全积分单元算法，其默认为单点积分算法。该单元可采用Lagrange列式，也可采用ALE列式，当采用ALE列式时，在关键字中需要定义单元算法关键字*SECTION_SOLID_ALE。单元支持大部分的LS-DYNA材料算法。BEAM161单元，为二维梁单元，用三个节点定义，I和J确定梁的轴向，K确定横截面主轴方位。具有Hughes-Liu及B-S算法，其默认算法为Hughes-Liu算法。单元可以定义不同的积分算法，默认值为22Guass积分。对于Hughes-Liu算法选项，如果截面积分规则Int.rule-arbitry.section取默认值0，则程序会采用梁单元跨中的一组积分点来模拟矩形及圆形的截面。 7.2.5材料失效方式的控制在有限元数值模拟中，由于程序本身是基于连续介质力学理论基础上形成的，因此，材料的失效是靠删除单元来实现的。单元达到屈服强度后，就被从模型中删除，不再参与计算，所以，对于材料失效强度的控制决定着数值模拟计算的准确性。材料失效现象比较复杂，因此，出现了各种失效理论模型，目前，比较成熟的强度理论模型主要有以下四种：第一强度理论模型，即最大拉应力理论模型； （3-1）第二强度理论模型，即最大伸长线应变理论模型； （3-2）第三强度理论模型，即最大剪应力理论模型； （3-3）第四强度理论模型，即形状改变比能理论模型（Von Mise应力）。 （3-4）四式中，-许用应力，=屈服极限/安全系数；，-三个主应力。对于铸铁、石料、混凝土和玻璃等脆性材料，通常以断裂形式失效，可采用第一和第二强度理论；对于碳钢、铜、铝等塑性材料，通常以屈服形式失效，可以采用第三和第四强度理论。根据以上论叙，对于混凝土材料，论文中采用材料控制项中FS和额外控制项*MAT_ADD_EROSION联合控制材料的失效。对于钢筋材料则采用材料模型*MAT_PLASTIC_KINEMATIC中提供的FS参数控制。7.2.6 接触方式的选择LS-DYNA提供丰富的接触方式可供选择，可以求解各种柔性体之间、刚性体之间以及柔性体和刚性体之间的接触问题，并可分析接触表面的静动摩擦力、固连失效以及流体与固体的界面等问题。程序中处理动态的接触-碰撞界面主要采用如下三种算法。（1）动力约束法动力约束法是最早采用的接触算法，该方法在每一步修正结构的构形之前，判断每一个没有与主表面接触的从节点是否会在该时步内贯穿主表面。如有从节点贯穿主表面，则将时步缩小，使那些从节点都不贯穿主表面，而其中有的刚到达主表面，在下一时步开始，对刚到达主表面的从节点施加动力约束条件。对所有已经与主表面接触的从节点都施加约束条件，以保持从节点与主表面的接触。此外，检查与主表面接触的从节点所属的单元是否存在受拉交界面力。如果受拉面力，则用释放条件使从节点脱离主表面。由于此算法比较复杂，后来只用于固连接触问题，即没有释放条件。（2）分配参数法该方法仅用于滑动界面接触问题的处理，其原理是：将所有正在接触的从单元质量的一半分配到主表面上，同时由每个单元的应力确定作用在接受质量的主表面面积上的分布压力。在完成质量和压力的分配后，程序修正主表面的加速度。然后对从节点的加速度和速度施加约束，以保证从节点沿主表面运动。这一算法不允许从节点穿透主表面，适合相对滑动但没有分离的表面之间的接触分析。（3）对称罚函数法对称罚函数法为LS-DYNA程序的缺省接触算法，该方法的基本原理为：每一时步先检查各从节点是否穿透主表面，没有穿透则对该节点不作任何处理。如果穿透，则在该从节点与被穿透主表面之间引入一个较大的界面接触力，其大小与穿透深度、主片刚度成正比，称为罚函数值。它的物理意义相当于在从节点和被穿透的主表面之间放置一个法向弹簧，以限制从节点对主表面的穿透。所谓对称罚函数法是指程序在对所有的主节点按如上步骤处理，其算法与从节点一样。对称罚函数法，很少激起网格的零能模式，没有噪声，这是由于算法具有对称性、动量守恒准确，不需要碰撞和释放条件。罚函数值大小受到稳定性限制。若计算中发现明显穿透，可放大罚函数值或缩小时间步长来调节。在LS-DYNA程序中，接触类型可以分为三类，即单面接触、节点-表面接触和表面-表面接触。（1）单面接触单面接触包括的接触类型有SS、ASSC、AG、ASS2D、ESS。单面接触可以用于一个物体表面各部分的自相接触或与另一个物体的表面接触，在单面接触中，程序会自动判定模型中发生表面接触的位置。因此，单元接触的定义是最简单的，无需定义源面和目标面，当定义好单面接触时，它允许一个模型的所有外表面都可能接触，有利于处理预先不知道接触表面的自身接触或大变形问题。（2）点-面接触LS-DYNA程序中的点-面接触类型有NTS、ANTS、RNTR、TDNS、TNTS、ENTS、DRAWBEAD、FNTS等。对于点-面接触，接触节点将穿透目标面，这种接触类型用于一般情况下的两个物体表面之间的接触。点-面接触需要定义target面和contact面，原则如下：平面或凹面为target面，凸面为contact面；网格较粗的面作为target面，网格较细为contact面。（3）面-面接触LS-DYNA程序中的面-面接触类型包括STS、OSTS、ASTS、ROTR、TDSS、TSTS、ESTS、SE、FSTS、FOSS、TSES等。当一个物体的表面穿透另一个物体表面时需用面-面接触，这种接触是最常用的接触类型，常用于任意形状且存在较大接触面积的物体之间的接触问题。该接触类型对于物体间有大量相对滑移时很有效。接触碰撞问题属于最困难的非线性问题之一，因为在接触碰撞问题中的响应是不平滑的。当发生碰撞时，垂直于接触界面的速度是瞬时不连续的。由于分析本身非常复杂，接触情况变化多端。因此，计算中选用了LS-DYNA提供的CONTACT-ERODING-SINGLE-FACE接触计算模型。该接触计算模型可以自动搜索接触面，判断接触，并可以处理侵蚀、断裂等复杂边界变化情况。接触刚度选用1.1，材料的摩擦系数统一设定为0.6。7.2.7重要参数的设定对于梁结构的每个端面和支撑柱的地面部分以及模拟地面的地板添加全方位的面约束。在生成关键字后，在关键字中改为Z方向的，重力加速度为9.8m/s。计算0到5秒内的爆破过程，将DATABASE_BINARY_D3PLOT设置为0.025s,将DATABASE_BINARY_D3THDT也设置为0.025s。由于计算机的计算能力限制每秒钟不能超过1000个plot文件，也就是说时间不能低于0.005s。7.2.8钢筋混凝土支撑梁结构爆破拆除数值模拟研究流程图钢筋混凝土支撑梁结构周围环境合理简化后建立实体模型选择合适的材料模型边界条件破坏准则选择合适的单元类型和单元确定网格尺寸并划分网格载荷施加设置计算时间生成关键字文件并修改计算图7.1钢筋混凝土支撑梁结构爆破拆除数值模拟研究流程图图7.1为本次数值模拟研究的流程图，整个建模以及计算过程都是在ANSYS12.0下完成的，建立的数值计算模型如图7.2所示。7.3.有限元模拟计算结果以及分析7.3.1爆破过程分析为实现对爆破设计中的延时起爆顺序的模拟，将支撑梁的不同爆区建立不同的部件，每个部件通过*MAT_ADD_EROSION来控制混凝土失效的时间从而来模拟爆破过程。每个爆区的延时时间为25ms,整个延时起爆过程如图7.3图7.20所示， 7.3.2混凝土单元的最大主应力随时间的变化此单元位于第一部分内，所以一开始就有一个压力突变，此单元为掉落在地上，是与钢筋连在一起，所以一直会有压力的变化。在0.4s时达到最大值。图表4-1此单元为第二部分的，在一开始会有强烈的压力突变，因为在中间部分，第二峰值是因为此单元掉落在地上，之后的波动则是单元之间的碰撞。第三部分，此单元紧挨立柱，所以一开始并无突变，在0.7s左右达到最大值，是因为单元一部分位于钢筋分开以至于掉落碰撞到住上，之后平稳。图表4-19图表4-18此单元在第四部分内，而属于中间部分，所以一开始就有一个小的压力突变在0.7s左右有强烈的突变是因为此单元掉落地上，之后的突变是因为她压在一个单元上面。第五部分内，此单元最后落在在支柱上面，且距离支柱的顶部有很短的距离，在一开始之后又两个突变的峰值，第一个是删除后的变形，第二是碰撞到支柱顶部的，之后不稳定是因为此单元还是有一部分连接着钢筋，切支柱是多个梁的交点。第六部分内，此单元与前边的分布距离较近，所以有两个相邻的峰值一个是前面部分删除引起的，一个是本身引起的。在后在0.8s左右的突变是掉落。图表4-21第七部分内，峰值是因为变形而产生的，没有一开始发生突变是因为此部分位于边沿，没有直接接触，且距离支柱较近，此面不是碰击地面的那部分，并且有一部分连接到钢筋上。第八部分内，此面在较长的一块混凝土上，与前一部分距离较近，所以有那两个压力突变一开始，在0.8s左右的压力突变是因为落地。有一部分在钢筋上。第九部分前面的突变是由之前的部分和自身部分共同作用引起的，之后的峰值是掉落在地上。后面的不平稳是因为虽然掉落在地上，但是还是有一部分在钢筋上。爆区II：第十部分，前面连续的变动是因为第一爆区，之后的突变是因为自身部分引起的，由于此梁整体较短，所以突变压力较小。之后的因为碰撞，未落地，有一部分与钢筋相连。图表4-25图表4-26第十一部分，之前突变与上面同理，之后的突变因为与地面碰撞。第十二部分，单元直接落地，后面两次压力突变是因为混凝土弹跳起来，两次碰撞地面。图表4-27第十四部分，同理上面因为弹起而两次触地，只是单元所在梁比之前的长一些，才会有不同。图表4-29第十六部分，此面与地面接触，最后的峰值是与地面接触造成的，此单元仍然与钢筋有连接部分。前面的压力突变是因为之前部分引起的。后面的波动是因为钢筋的抖动而使次面的压力发生变化。7.3.3钢筋单元的剪应力随时间的变化在第一部分内，且距离墙面的约束比较近，一开始就有突变，会发生剧烈的上下摆动，所以剪应力正负交替，最后在还是负的范围内波动。图表4-30在第二部分内，是支撑梁上部分的钢筋，在0.8秒达到峰值，发生剧烈变形。而且处于钢筋较密集的地方所以应力比较大。图表4-31在第三部分内，在钢筋和剩余混凝土的共同作用下，剪应力交替变化，变化幅度越来越小。距离支柱较远。图表4-33同样在第三部分内，钢此单元最后是由于混凝土压在上面，所以最后是负的剪应力，且变化不大。在0.5秒左右达到谷值，是因为混凝土下落接触，峰值是因为混凝土的弹起以及周边混凝土的作用。在第四部分内，属于上边的钢筋，所以在混凝土掉落时，没有直接的接触。图表4-35在第五部分内，在钢筋交叉部分附近，且附近有混凝土连接在钢筋上，在0.8秒左右达到峰值，是混凝土掉落而产生的力。在第六部分内的，在1秒左右有剧烈突变，是混凝土与钢筋的作用是此单元产生剧烈的变形。在第七部分内，是在混凝土中的