搜索
1混凝土力学性能数值仿真研究进展吉久茂，焦楚杰（广州大学土木工程学院，广东广州510006）摘要：回顾了目前研究混凝土数值仿真的常用方法，分别是有限元法，离散元法，边界元法及它们的发展过程。综述了数值仿真方法的基本原理和应用模型。侧重叙述了颗粒离散元法和散体细观力学理论及其DEM模型与改进。介绍了各种算法国内外研究进展及耦合算法的扩展研究情况。最后指出了一些数值仿真算法的存在问题，并对数值仿真的未来发展即科研实验与数值仿真结合研究提出了希望。关键词：混凝土数值仿真有限元离散元边界元1．数值仿真的起源与科研意义数值仿真技术也叫计算机模拟。它是以电子计算机为手段通过数值计算和图像显示的方法，达到对工程问题及物理问题研究的目的。数值模拟技术最早诞生于20世纪中叶。由国外知名学者BruceG.H和PeaceD.W模拟了一维气相不稳定径向和线形流。1955年，Peaceman与Rachford研发的交替隐式解法（ADI）取得了数值仿真技术的重大突破。该算法稳定且计算速度较快，所以在相关领域得到广泛应用。1959年，Peaceman与Douglas第一次进行了两维两相模拟，这标志着现代数值模拟技术的开始。60年代，CoatsK.H和NielsenR.L首次进行了三维两相模拟。1968年BreitenbachE.A发表了三维三相模拟解法，这期间的另一项成就就是Peaceman提出了后来通用的Peaceman方程，用于二维三维扩散方程的数值解答。80年代最大的成就是AppleyyardJR和CheshireI.M发表了嵌套因式分解法。90年代，ZoltanE.Heinemann提出了PEBI网格，PEBI网格结合了正交网格和角点网格的优点，现在正逐渐成为主流数值模拟网格体系。关于数值模拟，我们可以理解为利用计算机做实验，特别是近几十年来随着计算机性能的不断提高和力学学科的深入发展，数值模拟方法对问题的分析研究取得重大进展。特别是对数学分析方法难于处理的边界问题或复杂几何造型问题，以及对模拟试验方法代价过高的问题，数值模拟被证明是一种有效的解决问题的手段，为改善试验设计提供力学基础。————————国家自然科学基金项目（51278135，50708022）；广东省高等学校高层次人才项目（2050205）；建设部科研开发项目（06-K1-37、07-K4-5、07-K4-13、10-K3-27、10-K4-18）；广州市属高校“羊城学者”科研项目（10A043G）。第一作者简介：吉久茂（1989-），男，江苏南通人，硕士研究生，从事高强高性能混凝土问题研究。通讯作者简介：焦楚杰（1974-），男，湖南浏阳人，博士，教授，硕导，从事高强与高性能混凝土研究。（E-mail：jiaochujie@sina.com）22.有限元法的创立和发展概况在近代工程科学技术的发展中，由于飞行器，船舶，机械，水坝，桥梁，房屋等工程设计上的需要，固体力学始终受到人们的重视。20世纪40年代前虽然已提出变分法，差分法，松弛法等计算方法，但他们只能用于分析形状简单的结构，对于实际工程中很复杂的结构，事实上很难进行比较精确的分析。到了40年代中期，电子计算机出现以后，人们首先想到用计算机求解杆系结构力学中的力法和变位法的基本方程。形成矩阵力法和矩阵位移法。50年代中期，人们提出有限元法，把连续介质离散成一组单元，使无限自由度问题转化为有限自由度问题，再利用计算机进行求解。这种方法可以分析形状十分复杂的结构，所以受到专家学者的普遍重视，很快扩展到固体力学的各个分支。又从固体力学扩展到流体力学、热传导学、电磁学等各个领域，发展成为一个十分重要的工程计算方法。2.1有限单元法的基本原理有限单元法是一种获得工程问题近似解的数值分析方法，它通过离散的有限个单元集合体来代替真实的结构，通过对单元的分析，进而得到整个结构的近似解。先对预研究的结构进行单元划分，计算出各单元的单元刚度方程，然后组装各单元刚度的方程得到单元总刚。最后处理边界条件并求解得出节点位移，以及单元的应变应力等。在一般的工程设计中，这种近似解已经能满足实际工程的需要，有了足够的精度。有限单元法的最初形态是结构力学中所学的杆系结构矩阵位移法的一种推广，并应用于二维和三维的问题。然而与杆系结构不同，二维和三维实体没有明显的联接点。还需要建立许多人为的节点，将连续体离散化为许多任意形状的单元，用此方法，连续体便可用有限个自由度体系来近似表示，并求得其近似解。从求解的未知量来看，有限元法可分为应力场，位移场及混合场三种场变量。离散化的主要目标将物体分成充分小的单元，使得简单位移模型能足够近似精确解。在这个过程中，将人为通过网格划分出有限个单元。单元与单元之间以节点相连。作为位移场中有限单元的建立，我们所要求的是连接有限个单元体的节点的位移与节点力之间的转化关系，即确定单元的刚度矩阵。有限元矩阵的建立一般分为广义坐标有限模型和等参有限元模型。在实际分析中使用等参有限元模型更为有效，可以避免在广义矩阵中出现奇异矩阵的可能性，也减少了矩阵运算的次数。另外对曲边单元也能进行较为精确的模拟操作。然而广义坐标有限元模型的建立能更好的让我们加深对有限元法的理解。广义坐标法一般分为三类。即直接法，变分法，加权余数法。直接法易于理解，用于解决简单的问题。变分法把有限元归结为求泛函的极值问题。加权余数法直接从基本的微分方程出发，求出近似解。对于不存在泛函的工程领域，提供了有效的解决方法。2.2有限元法的优点有限元的优点是：①可以分析形状十分复杂、非匀质的各种实际问题的工程结构。②可以在计算中模拟各种复杂的材料本构关系、荷载和条件③可以对结构进行动力分析④由于前处理和后处理技术的发展可以进行大量方案的比较分析，并通过图形表示计算结果从而有利于对工程方案进行优化。3.离散元的起源及国内外的发展情况离散单元法的思想源于较早的分子动力学，是美国学者Cundall在七十年代率先提出的，并且广泛应用于岩石力学，土力学，结构分析等领域，是一种新兴的非连续体分析方法。1979年Cundall和strack提出适于土力学的离散元法，并推出二维圆盘程序BALL和三维圆球程序TRUBAL。1980年walton用此模型研究散体流动并有所进展。同年campbell提出硬颗粒模型并用于分析剪切流。1989年英国Aston大学Thornton引入了TRUBAC程序，能模拟干——湿，弹—塑性和颗粒两相流问题。澳大利亚新南威尔士大学余艾冰研究中心进行了多方面DEM模拟。CSIROC(commonwealthscientificandindustrialresearchorganization)研究所的3cleary用离散元模拟了不少工程问题。我国离散元研究始于20世纪80年代由王泳嘉引入学者cundall的离散元法进行岩石力学和颗粒系统的模拟。此后，块体元主要应用于边坡危害和矿井稳定等岩石力学问题当中。同济大学王卓琳博士基于离散单元法和刚体—弹簧法考虑混凝土由水泥砂浆，粗骨料以及二者界面层组成建立了混凝土的二维细观力学模型。通过单元间连接弹簧的破坏，表征混凝土材料开裂后由连续体向非连续体的转变。弹簧破坏后，单元间连接关系转化为接触关系。北京大学刘凯欣等利用离散元法，对冲击荷载下的混凝土平板的变形破坏全过程进行了数值模拟和动画显示，并分析了不同强制位移速率下不同的破坏形式。4.离散元法的单元类型离散元法的单元从形状上可以分为块体元和颗粒元两大类。块体元中最常用的有四面体元，六面体元。对于二维问题，还可以是任意多边形元。但应用范围不广，每个离散单元只有一个基本节点。也就是取其形心点。颗粒元主要是采用球体元，对于二维问题采用圆盘形单元。5.颗粒离散元法和散体细观力学理论散元法把离散体看作有限个离散单元的组合。根据其几何特征分为颗粒和块体两大系统。每个颗粒和块体为一个单元。根据过程中每一时步颗粒间的相互作用和牛顿运动定律的交替迭代预测散体群的行为。5.1运动学计算原理离散元法中颗粒运动学计算主要由在时步△t内球体颗粒的线运动和转动方程给定。Fi-ßgVi=m△Vi/△tMi-ßgWi=I△Wi/△t(i=1,2,3…)(1)其中Fi和Mi分别为不平衡力和力矩，Vi和△Vi为线速度及增量。Wi和△Wi为角速度及增量。m为质量，I为转动惯量。ßg为整体阻尼，解之可得速度，进而可得该球体新位置。Xi=Xi+Vi△t,Φi=Φi+Wi△t(2)由各球体的新位置坐标可决定相邻颗粒是否接触或原接触体是否分离。相互接触的球体会产生假性重叠，再由接触模型公式分别求出接触力Fi和Mi，返归（1）式进行迭代。5.2接触模型接触模型是颗粒离散元法的核心。颗粒离散元中主要分两种模型。一种是干颗粒模型，另一4种是湿颗粒模型。干颗粒模型是接触的两圆球单元在法-切向相对运动时接触力和局部变形的拟静态关系。有许多研究者仍采用弹簧-阻尼器模型，根据经验或实验给定参数。湿颗粒模型是一种近似模型。当两球形单元对心相对运动时，流体黏性会产生挤压力以及切向相对运动时的阻力。两种模型的接触力和变形关系都是非线性的。法-切向作用很难分开，但是单一的法向或切向求解又很复杂，故近似采用叠加原理。5.3.DEM颗粒接触模型Cundall首创的盘元接触模型如下图。图中Kn为法向刚度，Kt为切向刚度。dn为法向阻尼，dt为切向阻尼。因为该模型相对简单，故现今仍有许多人采用，也有许多专家学者在此基础上进行了改进。（a）法向作用(b)切向作用（c）滚转作用5.3.1DEM改进模型有些学者在cundallDEM模型的基础上提出若干改进接触模型。（1）Oda的改进离散元法，考虑法向接触应力不对称性，简化出接触力矩即在cundall的提出模型的基础上不仅考虑了法向和切向作用还考虑了单元间的滚动作用。（c）图中Kr为滚转刚度，dr为滚转阻尼。（2）Kishino的粒状元素法。基本模型和cundall相同。主要特点是用接触点构成的瞬时刚度矩阵控制单元运动。适合于应力控制下的拟静态计算。5.3.2DEM其他计算模型国内专家的研究进展北京大学力学与工程科学系刘凯欣等提出了三维连结型离散元模型，建立了可实现连结型—接触型转化的三维离散元计算程序，用于模拟由连续介质转化为非连续介质的过程。该模型是将被解空间离散成按一定的空间位置排列的刚性球体离散单元的结合，单元间由一个法向和两个切向方向的三种线型弹簧相互连接。弹簧的刚度根据连续介质力学的基本原理得到。每个单元中心位置的变化是受到单元间弹簧的作用力所控制。在此基础上，建立了三维离散元计算程序，用来模拟连续介质转变为非连续介质的力学过程。利用该模型仿真了冲击荷载下材料应力波传播过程。将计算结果与LS-DYNA计算结果进行了比较，结论大致相同。表明了离散元模型不但可以应用于连续介质动力响应问题，也可应用于非连续介质力学问题。浙江大学金伟良等提出了矩形离散单元法模型。通过引入节点单元和考虑混凝土的非线性对矩形离散单元进行了进一步改进，并给出了改进后离散元模型的破坏的准则和基本方5程以及弹簧系数，改进模型采用双链表技术，提高了模型的计算效率。同济大学吕西林等为了研究钢筋混凝土结构在地震作用荷载下的力学性能，建立了钢筋混凝土杆件的杆端多弹簧模型。杆端多弹簧模型将构件沿纵向分割成n个杆段，每一个杆段为一个不发生变形的刚体单元，刚体单元之间由连接刚体单元质心截面的轴向弹簧组相连。轴向弹簧组代表代表相邻质心长度范围内构件的力学性能。轴向弹簧组包含数个混凝土弹簧与钢筋弹簧。将构件截面沿横向分成m条带，混凝土弹簧位于条带面积中心，代表条带混凝土力学性能。钢筋弹簧位于条带内钢筋的中心位置，代表条带范围内钢筋的力学性能。6.离散元与其他数值方法耦合研究进展离散元法与其他数值方法结合进行扩展性研究是一个热点。如离散元与边界元结合研究岩石力学，与计算流体动力学结合研究颗粒两相流。与有限元结合研究结构工程问题。近年来，有限元法和离散元法结合解决工程问题的最新进展是Han等的喷丸成形模拟和Dwen等的多重断裂固体和离散系统的模