混凝土力学数值仿真研究进展 (修改)

混凝土力学性能数值仿真研究进展吉久茂 ，焦楚杰（广州大学 土木工程学院，广东 广州 510006）摘 要：回顾了目前研究混凝土数值仿真的常用方法，分别是有限元法，离散元法，边界元法及它们的发展过程。综述了数值仿真方法的基本原理和应用模型。侧重叙述了颗粒离散元法和散体细观力学理论及其DEM模型与改进。介绍了各种算法国内外研究进展及耦合算法的扩展研究情况。最后指出了一些数值仿真算法的存在问题，并对数值仿真的未来发展即科研实验与数值仿真结合研究提出了希望。关键词：混凝土 数值仿真 有限元 离散元 边界元71数值仿真的起源与科研意义数值仿真技术也叫计算机模拟。它是以电子计算机为手段通过数值计算和图像显示的方法，达到对工程问题及物理问题研究的目的。 数值模拟技术最早诞生于20世纪中叶。由国外知名学者Bruce G.H和PeaceD.W模拟了一维气相不稳定径向和线形流。1955年，Peace man与Rach ford 研发的交替隐式解法（ADI）取得了数值仿真技术的重大突破。该算法稳定且计算速度较快，所以在相关领域得到广泛应用。1959年，Peace man与Douglas第一次进行了两维两相模拟，这标志着现代数值模拟技术的开始。60年代，Coats K.H和Nielsen R.L首次进行了三维两相模拟。1968年Breitenbach E.A发表了三维三相模拟解法，这期间的另一项成就就是Peace man提出了后来通用的Peace man方程，用于二维三维扩散方程的数值解答。80年代最大的成就是Appleyyard J R和Cheshire I.M发表了嵌套因式分解法。90年代，Zoltan E.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）2.有限元法的创立和发展概况在近代工程科学技术的发展中，由于飞行器，船舶，机械，水坝，桥梁，房屋等工程设计上的需要，固体力学始终受到人们的重视。20世纪40年代前虽然已提出变分法，差分法，松弛法等计算方法，但他们只能用于分析形状简单的结构，对于实际工程中很复杂的结构，事实上很难进行比较精确的分析。到了40年代中期，电子计算机出现以后，人们首先想到用计算机求解杆系结构力学中的力法和变位法的基本方程。形成矩阵力法和矩阵位移法。50年代中期，人们提出有限元法，把连续介质离散成一组单元，使无限自由度问题转化为有限自由度问题，再利用计算机进行求解。这种方法可以分析形状十分复杂的结构，所以受到专家学者的普遍重视，很快扩展到固体力学的各个分支。又从固体力学扩展到流体力学、热传导学、电磁学等各个领域，发展成为一个十分重要的工程计算方法。2.1有限单元法的基本原理有限单元法是一种获得工程问题近似解的数值分析方法，它通过离散的有限个单元集合体来代替真实的结构，通过对单元的分析，进而得到整个结构的近似解。先对预研究的结构进行单元划分，计算出各单元的单元刚度方程，然后组装各单元刚度的方程得到单元总刚。最后处理边界条件并求解得出节点位移，以及单元的应变应力等。在一般的工程设计中，这种近似解已经能满足实际工程的需要，有了足够的精度。有限单元法的最初形态是结构力学中所学的杆系结构矩阵位移法的一种推广，并应用于二维和三维的问题。然而与杆系结构不同，二维和三维实体没有明显的联接点。还需要建立许多人为的节点，将连续体离散化为许多任意形状的单元，用此方法，连续体便可用有限个自由度体系来近似表示，并求得其近似解。 从求解的未知量来看，有限元法可分为应力场，位移场及混合场三种场变量。离散化的主要目标将物体分成充分小的单元，使得简单位移模型能足够近似精确解。在这个过程中，将人为通过网格划分出有限个单元。单元与单元之间以节点相连。 作为位移场中有限单元的建立，我们所要求的是连接有限个单元体的节点的位移与节点力之间的转化关系，即确定单元的刚度矩阵。有限元矩阵的建立一般分为广义坐标有限模型和等参有限元模型。在实际分析中使用等参有限元模型更为有效，可以避免在广义矩阵中出现奇异矩阵的可能性，也减少了矩阵运算的次数。另外对曲边单元也能进行较为精确的模拟操作。然而广义坐标有限元模型的建立能更好的让我们加深对有限元法的理解。广义坐标法一般分为三类。即直接法，变分法，加权余数法。直接法易于理解，用于解决简单的问题。变分法把有限元归结为求泛函的极值问题。加权余数法直接从基本的微分方程出发，求出近似解。对于不存在泛函的工程领域，提供了有效的解决方法。2.2有限元法的优点有限元的优点是：可以分析形状十分复杂、非匀质的各种实际问题的工程结构。可以在计算中模拟各种复杂的材料本构关系、荷载和条件可以对结构进行动力分析由于前处理和后处理技术的发展可以进行大量方案的比较分析，并通过图形表示计算结果从而有利于对工程方案进行优化。3.离散元的起源及国内外的发展情况 离散单元法的思想源于较早的分子动力学，是美国学者Cundall在七十年代率先提出的，并且广泛应用于岩石力学，土力学，结构分析等领域，是一种新兴的非连续体分析方法。1979年Cundall和strack提出适于土力学的离散元法，并推出二维圆盘程序BALL和三维圆球程序TRUBAL。1980年walton用此模型研究散体流动并有所进展。同年campbell提出硬颗粒模型并用于分析剪切流。1989年英国Aston大学Thornton引入了TRUBAC程序，能模拟干湿，弹塑性和颗粒两相流问题。澳大利亚新南威尔士大学余艾冰研究中心进行了多方面DEM模拟。CSIROC(common wealth scientific and industrial research organization)研究所的cleary用离散元模拟了不少工程问题。我国离散元研究始于20世纪80年代由王泳嘉引入学者cundall的离散元法进行岩石力学和颗粒系统的模拟。此后，块体元主要应用于边坡危害和矿井稳定等岩石力学问题当中。同济大学王卓琳博士基于离散单元法和刚体弹簧法考虑混凝土由水泥砂浆，粗骨料以及二者界面层组成建立了混凝土的二维细观力学模型。通过单元间连接弹簧的破坏，表征混凝土材料开裂后由连续体向非连续体的转变。弹簧破坏后，单元间连接关系转化为接触关系。北京大学刘凯欣等利用离散元法，对冲击荷载下的混凝土平板的变形破坏全过程进行了数值模拟和动画显示，并分析了不同强制位移速率下不同的破坏形式。4.离散元法的单元类型离散元法的单元从形状上可以分为块体元和颗粒元两大类。块体元中最常用的有四面体元，六面体元。对于二维问题，还可以是任意多边形元。但应用范围不广，每个离散单元只有一个基本节点。也就是取其形心点。颗粒元主要是采用球体元，对于二维问题采用圆盘形单元。5.颗粒离散元法和散体细观力学理论 散元法把离散体看作有限个离散单元的组合。根据其几何特征分为颗粒和块体两大系统。每个颗粒和块体为一个单元。根据过程中每一时步颗粒间的相互作用和牛顿运动定律的交替迭代预测散体群的行为。5.1运动学计算原理离散元法中颗粒运动学计算主要由在时步t内球体颗粒的线运动和转动方程给定。Fi-gVi=mVi/tMi-gWi=IWi/t (i=1,2,3) (1)其中Fi和Mi分别为不平衡力和力矩，Vi和Vi为线速度及增量。Wi和Wi为角速度及增量。m为质量，I为转动惯量。g为整体阻尼，解之可得速度，进而可得该球体新位置。Xi=Xi+ Vit, i=i+ Wit (2)由各球体的新位置坐标可决定相邻颗粒是否接触或 原接 触体是否分离。相互接触的球体会产生假性重叠，再由接触模型公式分别求出接触力Fi和Mi，返归（1）式进行迭代。5.2接触模型接触模型是颗粒离散元法的核心。颗粒离散元中主要分两种模型。一种是干颗粒模型，另一种是湿颗粒模型。干颗粒模型是接触的两圆球单元在法-切向相对运动时接触力和局部变形的拟静态关系。有许多研究者仍采用弹簧-阻尼器模型，根据经验或实验给定参数。湿颗粒模型是一种近似模型。当两球形单元对心相对运动时，流体黏性会产生挤压力以及切向相对运动时的阻力。两种模型的接触力和变形关系都是非线性的。法-切向作用很难分开，但是单一的法向或切向求解又很复杂，故近似采用叠加原理。5.3.DEM颗粒接触模型Cundall首创的盘元接触模型如下图。图中Kn为法向刚度，Kt为切向刚度。d n 为法向阻尼，d t 为切向阻尼。因为该模型相对简单，故现今仍有许多人采用，也有许多专家学者在此基础上进行了改进。（a）法向作用 (b)切向作用 （c）滚转作用5.3.1 DEM改进模型有些学者在cundall DEM模型的基础上提出若干改进接触模型。（1） Oda的改进离散元法， 考虑法向接触应力不对称性，简化出接触力矩即在cundall的提出模型的基础上不仅考虑了法向和切向作用还考虑了单元间的滚动作用。（c）图中Kr为滚转刚度，d r为滚转阻尼。（2） Kishino的粒状元素法。基本模型和cundall相同。主要特点是用接触点构成的瞬时刚度矩阵控制单元运动。适合于应力控制下的拟静态计算。5.3.2 DEM其他计算模型国内专家的研究进展北京大学力学与工程科学系刘凯欣等提出了三维连结型离散元模型，建立了可实现连结型接触型转化的三维离散元计算程序，用于模拟由连续介质转化为非连续介质的过程。该模型是将被解空间离散成按一定的空间位置排列的刚性球体离散单元的结合，单元间由一个法向和两个切向方向的三种线型弹簧相互连接。弹簧的刚度根据连续介质力学的基本原理得到。每个单元中心位置的变化是受到单元间弹簧的作用力所控制。在此基础上，建立了三维离散元计算程序，用来模拟连续介质转变为非连续介质的力学过程。利用该模型仿真了冲击荷载下材料应力波传播过程。将计算结果与LS-DYNA计算结果进行了比较，结论大致相同。表明了离散元模型不但可以应用于连续介质动力响应问题，也可应用于非连续介质力学问题。浙江大学金伟良等提出了矩形离散单元法模型。通过引入节点单元和考虑混凝土的非线性对矩形离散单元进行了进一步改进，并给出了改进后离散元模型的破坏的准则和基本方程以及弹簧系数，改进模型采用双链表技术，提高了模型的计算效率。同济大学吕西林等为了研究钢筋混凝土结构在地震作用荷载下的力学性能，建立了钢筋混凝土杆件的杆端多弹簧模型。杆端多弹簧模型将构件沿纵向分割成n个杆段，每一个杆段为一个不发生变形的刚体单元，刚体单元之间由连接刚体单元质心截面的轴向弹簧组相连。轴向弹簧组代表代表相邻质心长度范围内构件的力学性能。轴向弹簧组包含数个混凝土弹簧与钢筋弹簧。将构件截面沿横向分成m条带，混凝土弹簧位于条带面积中心，代表条带混凝土力学性能。钢筋弹簧位于条带内钢筋的中心位置，代表条带范围内钢筋的力学性能。6.离散元与其他数值方法耦合研究进展 离散元法与其他数值方法结合进行扩展性研究是一个热点。如离散元与边界元结合研究岩石力学，与计算流体动力学结合研究颗粒两相流。与有限元结合研究结构工程问题。近年来，有限元法和离散元法结合解决工程问题的最新进展是Han等的喷丸成形模拟和Dwen等的多重断裂固体和离散系统的模拟研究。能预测工程中连续体在一定条件下部分破裂-破碎以及爆炸力学问题。7.边界元法的基本思想在数值模拟处理边界问题的工作中往往采用多种计算方法耦合的算法。边界元也是新兴的一种数值仿真处理工具。物理问题的数学模型的形成，一般有几种不同的形式，可以直接表示成微分方程的形式，也可表现为区域上的变分形式。也可归结为边界方程上某个积分方程的形式。这些不同的数学形式在理论上是等价的。边界元方法也称为边界积分方程法. 它的基本思想是基于格林公式的应用一个区域上,积分转化为区域边界上的积分。早期的边界积分方程法, 只是指把偏微分方程的边值问题归化为边界积分方程的公式化过程。近代的边界积分方法, 不仅包含了各种方式的边界归化过程, 更重要的是包含了求边界积分方程数值解的各种形式的离散化技术。英国学者Brebbia提出作为一种与有限差分法，有限单元法并列的离散化数值方法，边界积分方程法应该作为边界元法成为相对独立的一个数学方法分支。从此边界元的说法便在学术界流传开来。边界元法是在经典的边界积分方程法的基础上吸收了有限元离散化技术而发展起来的一种偏微分方程的数值解法. 它的主要优点是可将空间的维数降低一维, 且只要求出边界上定义在所谓节点上的一定的值(虚拟密度, 位势和流量, 位移和面力等)另外弹性力学中的平衡微分方程，几何方程，应变协调方程，广义胡克定律，边界条件等方程的建立也为有限元、边界元的数值计算提供了基础。8.边界元法耦合算法在数值模拟方法中往往也有耦合算法。用耦合法的本质其实就是为了扬长避短。把一个任意的求解区域分成两个或多个子区域。问题可以使非线性非均匀性，尽管这个子区域不是相对规则的，但运用有限元却是方便的；在另一个含奇异性可以是无限的但是在几何上是规则的子区域上，问题转变为线性的均匀的，正好在这个子区域上采用边界元方法。两者相结合，就可以得到整个区域上的解。其耦合方法还可进一步用于解决无界区域上的平面弹性问题及stokes方程组等。从数值计算的角度来说，比较常用的是将使其易于与现有的有限元程序相结合。但应当防止总刚度矩阵的对称性破坏及带宽的增加等问题。另一种耦合的途径，是由有限元方程产生一系列边界条件用于边界元的公式中，以实现有限元与边界元的结合得到边界元类型的线性方程组。边界元与离散元耦合算法国内也有相关专家研究。清华大学金峰等选择动力时域边界元法与可变形体离散元耦合，提出了二维时域动力边界元-可变形体离散元耦合模型，充分发挥离散元与边界元的优点，将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中，为地下结构和岩质边坡的抗震稳定分析提供了全新的手段。9.边界元法的不足边界元方法也有其弱点。它对变系数非线性问题的适应性不如有限元法。当所考虑的区域的表面积与其体积之比值较大时, 边界元法不够有效。 在数值计算方面, 由于积分核的奇异性和离散化后的线性代数方程组的系数矩阵的非稀疏性和非对称性也增加了一些困难。10.结束语我们了解到有限元法，离散元法，边界元法等数值模拟方法既有他们独特点，也存在不足之处。如果我们将理论分析，试验方法，数值方法相结合，互为补充，必将取得相得映彰的效果。 参考文献：（1）龙驭球，有限元法概论M 北京: 人民教育出版社（2） 王强,吕西林. 离散单元法及其在建筑工程中的应用现状A 。 现代土木工程理论与实践C 。南京:河海大学出版社,2003. 656 660.（3）刘凯欣，高凌天，郑文刚. 混凝土动态破坏过程的数值模拟。 工程力学, 2000:470474.（4）祝家麟 边界元方法一个老而新的发展中的数值方法 重庆建筑工程学院 知识与进展（5）刘凯欣，高凌天 离散元法研究的评述, 力学进展 （6）金伟良，方韬 钢筋混凝土框架结构破坏性能的离散单元法模拟 工程力学 2005 第22卷 第4期（7）刘凯欣，高凌天，郑文刚 离散元法在求解三维冲击动力学问题中的应用 固体力学会议论文集（8）金峰，贾伟伟，王光纶 离散元-边界元动力耦合模型 水利学报 2001年（9）王强,吕西林，雷淑忠 离散单元法在钢筋混凝土构件非线性分析的应用 沈阳建筑大学学报 自然科学版（10）方韬，龚顺风，金伟良 混凝土结构破坏过程的离散单元法模拟 浙江大学学报 （工学版）（11）Day A S.A Introduction to dynamic relaxtion J.The engineer 219（3）：218224（12）Scott B D,Parck R,Priestley M J N.Stress-strain behavior of concrete confined overlapping hoops at low and high strain ratesJ.Microcomputers in civil Engineering 1992,7（2）：152159（13）Iwashita K Oda M.Micro-deformation mechanism of shear banding process based on modified distint element method. Powder Technology ,2000,109:192205（14）Tsuji Y.Activities in discrete partical simulation in Japan. Powder Technology ,2000,113（15）Yang S C,Hsiau S S.Self-diffusion analysis in a vibrated granular bed Advanced Powder Technology 2001（16）Cundall P A A computer model for simulating prograsive large scale movement in block rock systerm.Symposium ISRM,1971,Proc 2（17）Tanaka H,Momozu M,Oida A,Yamazaki M.Simulation of soil deformation and resistence at bar penetration by the distinct element method. Journal of Terramechnics,2000（18）Cleary ,P W .predicting charge motion,power draw ,segregation and wear in ball mills using discrete element method.Nuclear Engeneering,1998（19）Mansur WJ ;Brebbia CA Topics in Boundary Element Research 87123,1985（20）Brebbia C.A,The Boundary Element Method for engineers,Pentech Press ,London,1981（21）朱伯芳 有限单元法原理与应用 中国水利水电出版社 知识产权出版社（22）王泳嘉，邢纪波 离散元法及其在岩土力学中的应用 辽宁 东北工学院出版社 1991 6089（23）邢纪波，王泳嘉 离散元法的改进及其在颗粒介质研究中的应用 岩土工程学报 1990 12（5）；5157（24）焦楚杰，蒋国平 钢纤维高强混凝土抗冲击性能的数值仿真 江苏大学学报（自然科学版） 2012年5月 第3卷 第33期 351353（25）Han k,peric D,Owen D R J,Yu J.A Combined finite /dicrete element simulation of shot peening processes part;3d interation laws Engng Comp,2000,17(6-7) 680702（26）Hank,Owen D R J,Peric D Combined finite/discrete element and explicit/implicit simulations of peeb forming processes.Engng Comp,2002,19(1);92118（27）Owen D R J,Feng Y T.Parallelised Finite/Discrete Element simulation of Mult-fracture solids and Discrete systems. Engng Comp,2001,18(3-4);557576Research progress of numerical simulation of the mechanical properties of concreteJI Jiumao JIAO Chujie(School of Civil Engineering，Guangzhou university，Guangzhou 510006，China)Abstract: Review of current research used concrete numerical simulation method, the finite element method, discrete element method, the boundary element method and their development process. The basic principles and applications of the model for the numerical simulation method. Focused narrative micromechanics particle discrete element method, and the dispersion of theory and its DEM model. Extension study of various algorithms research advances and coupling algorithm. Concluded that the existence of the problem of the numerical simulation algorithm, and the numerical simulation of the future development of scientific experiments and numerical simulation study hope.Key Words: discrete element boundary element the concrete numerical simulation finite element