大坝混凝土三维细观力学数值模型研究pdf

第25卷第2期计算力学学报V01．25，No．22008年4月ChineseJournalofComputationalMechanicsApril2008文章编号：1007—4708(2008)02一0241一07大坝混凝土三维细观力学数值模型研究马怀发“，陈厚群1，吴建平2，黎保琨3(1．中国水利水电科学研究院工程抗震中心，北京100044}2．并行与分布处理国家重点实验室，长沙410073；3．北京工业大学，北京100044)摘要：在细观结构层次上将大坝混凝土作为骨料、固化水泥砂浆及其粘结界面组成的复合材料．建立了大坝混凝土三维细观力学数值模型。该模型既能够反映混凝土及其细观各相材料在荷载作用下的损伤演化过程．又考虑了动载作用的应变率强化效应．给出了该数值模型求解方法，并编制出能够在普通PC机上运行的串行程序。加载过程既可按荷载控制又可按位移控制。同时，为了减少求解自由度应用了分尺度方法以使最小骨料和固化水泥砂浆混合后其力学性能与一种复合介质等效。通过混凝土湿筛和三级配试件的静、动(冲击)弯拉数值计算验证了本文计算方法和程序正确有效。另外．在串行程序的基础上，优化了刚度矩阵的存储方式，采用双门槛不完全Cholesky分解(ICT)预条件的共轭梯度法(CG)，完成了能够在SunFire6800服务器实现并行计算的并行程序改造。从而大大提高了计算效率．关键词：大坝混凝土；三维随机骨料模型，分尺度方法；预条件；并行计算中图分类号：TV313文献标识码：A．j．一J引置混凝土是由水、水泥和粗细骨料组成的复合材料。骨料尺寸及其级配、水泥用量、水灰比、固化强度、养护条件、环境湿度和混凝土的发热量等因素直接影响混凝土力学特性。不均匀性是混凝土材料的本质特点，微裂缝是决定其性能的主导因素。为了反映混凝土细观结构不均匀性特征，对其损伤断裂过程进行细观力学数值模拟，人们提出了许多细观力学模型，如格构模型，随机力学特性模型和随机骨料模型等。但到目前为止，对混凝土细观力学的数值模拟研究工作主要限于平面问题[1-]。尽管有三维骨料模型的研究[5]，但基本上停留在对“凸型”多面体骨料随机投放算法的研究，并且骨料含量低，没有考虑实际骨料级配。文献[6]提出了按照混凝土实际级配的骨料含量和比例生成球形骨料模型的随机投放方法，并在球形骨料模型的基础上生成了接近于实际骨料形态的多面体骨料模型。本文将在其基础上，针对混凝土三维细观力学收稿日期：Z006—03—15f修改稿收到日期：2006—10—16，基金项目：国家自然科学基金(90510017)重点资助项目．作者简介：马怀发’(1962·)，男，高级工程师(E-mail；mahf(园iwhr．tom)l陈厚群(1932一)。男．中国工程院院士．数值分析方法及其编程技巧进行分析研究。2基本方程2．1动力学方程混凝土及其细观各相材料的抗拉强度的强化系数H。和弹性模量强化系数H￡与应变率的关系表示为H，=exp{[A，(19㈦+B，)]c，)魄一exp{[A￡(19l三I+B￡]％))(1)式中A，，且，Cf为强度强化参数；AE，BE，CE为弹性模量强化参数。动拉强度、动弹性模量可表示为五(三)=H，Z，，E(三)=H￡E，(2)式中正，为混凝土材料的静态抗拉强度，E为静态弹性模量。^(；)及EG)表示混凝土动抗拉强度和动弹性模量，是拉(或压)应变率；的函数，并认为泊松比不受率效应的影响。动力学增量方程‘73为[明(△机(￡)}+[c](△抚(￡)}十[KJ(￡)](AUa(￡))=(△Pd(￡)}+{△P：(D))+{△P，(D)}(3)式中[明为质量阵，[C]为阻尼阵，[圮(￡)]表示万方数据计算力学学报第25卷动刚度阵；{△吼(￡))，(△讥(￡))，{AUa(t)}分别为结点加速度增量、速度增量和动位移增量，(△R(￡)}为结点动荷载列阵增量。{AP'd(D)}为在动位移上因损伤产生的应力转移而产生的附加荷载，{△P，(D))为在静位移上因损伤产生的应力转移而产生的附加荷载。采用文献[3]给出的双折线弹性损伤演化模型，残余强度系数A=^，／^，^为混凝土及其细观各相材料的抗拉强度，厶破坏单元的抗拉残余强度；与残余强度对应的应变系数叩一Cr／e。，e。为单元应力达到抗拉强度时的主拉应变，er为与抗拉残余强度相对应的最小应变；极限应变系数e=￡。／￡。，钆为极限拉应变。2．2方程求解步骤按照纽马克插值方法，取参数艿≥0．5，y=0．25(O．5+艿)2，令积分常数：】艿1ao。—：YA—t2’吼一涵’眈2涵口’2。历=历’1，铂口。2=7’争，(25=虿譬I(歹争一zJ2)(静力计算时，令口o=口l=口2=a3=口4=口5=0)将动力学方程(3)离散为[K。+aoM+口1C]△批(￡)=△PJ(￡)+AP'd(￡)+△只(￡)+M[a2讥+口3Ut]+C[口4以+口5ut]下面是求解该动力学离散方程求解步骤：(1)给出初值(2)选择时间步长、参数’，和艿，并计算积分常数(3)对咒+1时间步，即在t+At时间步内①初始化i=0K州=K。(D跪)形成刚度矩阵；Auo+。=0初始化位移增量；“o州=‰初始化总位移；矗：+，=矗。初始化速度，’“o。0+。=ii。初始化加速度②进行第i次迭代：(a)计算方程右边项：△P—APd。+l+AP；什l+△P。+l+^彳(口2矗。+a3瑟。)+C(口4矗。+n5砬。)△P0。一∑HE(k)AD卜1(￡)KoUd卜1△P科l=∑ADH(￡)K。∥式中Ko为初始刚度矩阵，损伤参数增量ADH(￡)一D(ei)一D(e。)(b)求位移增量：[K。+aoM+口lC]Au‘=AF‘Au‘=[K。+aoM+口1C]qAF‘(c)总动位移暂取为“算。=‰+Aui(d)计算每个单元的应力应变状态：按动拉强度判断是否有新的单元破坏，并计算新破坏单元的损伤参数；求应变增量△e0。=BAui，求动应变e11=e。+△￡0l，总应变￡井T(Dl=e’+￡井il，并由￡鬻求总主应变，由总主应变考察已受损单元是否有进一步的损伤，求进一步损伤单元的损伤参数D(e0。)，并求刚度增量矩阵(e)验证是否收敛：位移收敛准则为np∑0△彰一△扰芦1l|≤Pr。』摹1式中PrD为满足平衡要求的位移容许误差参数。若不满足精度要求，迭代次数小于最大容许次数进行下一步迭代，迭代次数大于最大容许次数，或若满足精度要求转入第(f)步(f)计算本荷载步的动应力％+。=(DoH辑。(钆z辑。，静应力Z+。一D2￡‘，动位移“椭=‰+Au‘，速度矗。+1=矗。+口lAu‘一口4h。一口5玩，加速度娩抖l=玩+aoAu‘一口2讥一口3玩，损伤参量D(￡州)=D(E。)+△D(￡：+1)．(g)总应力矗。=盯州+O-s和总位移“T+。=“州+“‘，并由矗。求总主应力。(h)根据主应变率计算弹性模量强化参数HE(三。+。)和抗拉强度强化参数H，(毒时。)．(4)判断试件是否产生失稳破坏或已完成给定时间步，输出结果计算停止，否则进行下一时间步计算。在程序中，采用了动态数组存储方式，模块化编程技巧，方程求解采用了克劳特(Crout)三角分解法。为了使程序能够在普通PC机上运行，在三维程序设计中采用了整体刚度矩阵分块求解算法，外存储读写方式。3三维随机骨料模型数值计算抗弯数值试验为三分点加载方式咖。采用随机球骨料模型。混凝土及其各相组分材料力学特性参万方数据第2期马怀发，等：大坝混凝土三维细观力学数值模型研究表1混凝土各相组分材料力学特性参数Tab．1Materialparametersofmeso—structure茁ofconcrete至＼材料1再i％泊松比1再i淼容重／kN·—n一3采簧黼蜮蝴舭∞拍n孔∞∞∞∞∞％帖∞加∞∞撕∞∞巧他嘶∞昨∞盯n孔孔∞∞∞∞损伤及强化参敷度系数变系数变系数强度强化参数弹模强化参数★’{AlB|c|AEBEcE0．103．0010．000．156．002．000．106．002．00数取值列入表1，认为细观各相单元的抗拉强度和弹性模量遵循对数正态分布[4]。在应力一位移曲线中，取混凝土试件下边缘中点的名义应力(仃一PL／6^2，P为荷载(／D为竖向位移)，L为支座间距，6为试件厚度，h为试件高度)与试件上边缘中点竖向位移(v／h)相对应。弯拉强度为应力一位移失稳时极限名义应力。在动力计算过程中，采用瑞利阻尼，阻尼比取0．02，抗拉强度、弹性模量的率效应分别按式(1)和式(2)计算。3．1湿筛试件细观力学分析试件尺寸为150mmx150mm×550mm，取跨中宽度为150mm为细观剖分区域，其余部分单元取宏观均质混凝土力学特性参数。豳Z湿筛混凝土弯拉细观效值试验Fig．2Theflexuralfailureofnumericalsimulationofmeso-modelofwet-siftedconcretespecimen骨料约为混凝土总重量的47％，其中，中骨料粒径为30mm，45颗；小骨料粒径为15mm，438颗。将细观区域剖分为6面体单元，其棱边长为3．75mm，为小骨料粒径的四分之一。用4面体单元将细网格6面体单元过渡到两边粗网格6面体单元。整个剖分区域共有76785节点，112020个单元，网格剖分如图1所示。计算得到应力一位移盟线如图2所示。静力分析有两种加载方式，按荷载控制时，静加载步长为0．083kN，应力一位移曲线为静力分析(1)：静弯拉强度为3．72MPa；按位移控制，竖向位移增量为0．001mm，应力一位移曲线为静力分析(2)：横轴仍为试件上边缘中点竖向位移(臼／^)，而纵轴为有限元计算得到的混凝土试件下边缘中点节点的．．717向正应力的均值。得到静弯拉强度为3．84MPa。动加载速率为168kN／s，按荷载控制加载，动弯拉强度约为4．66MPa。在文献[9]中，试验测得的静弯拉强度均值3．84MPa，动弯拉强度均值为4．73MPa。计算值接近于试验值。图3静极限荷载作用下混凝土试件弯折损伤云图图l湿筛混凝土试件细观模型Fig．3FlexuraldamagecontoursofconcretespecimenFig．1Meso-scalenumericalmodelofwet-siftedconcretespecimenwithultimatestaticloading万方数据计算力本文程序可以生成与有限元前后处理软件GID接口的文件。通过GID软件平台可以云纹图、等值线图、矢量图、变形图、曲线、动画等形式输出计算得到的位移场，应力场，单元损伤参数以及弯折损伤失稳过程。图3给出了静载弯折失稳时混凝土试件变形和损伤云图。3．2小骨料的“熔化”由于细观有限元网格划分中的单元尺寸一般应小于最小骨料粒径的1／3～1／4，尤其是界面层很薄，其单元尺度更小，这样节点自由度达到几十万甚至上百万。要进行大尺度全级配试件三维数值模拟，必须解决计算速度这个瓶颈问题。为了缩短计算时间，一种办法是使用大型高性能计算机并开发并行计算程序，另一种思路就是下面要讨论的方法，即利用分尺度理论[10]简化三维细观计算模型。试验结果表明，在大坝混凝土的弯拉试验中，破坏主要是沿最为薄弱的中、粗骨料和水泥砂浆的界面开始起裂，逐渐向水泥砂浆介质中发展。只有在等级较高的混凝土中，其骨料为人工破碎且质地较差时，裂缝有可能会延伸到中、粗骨料。但通常粒径在(5～20)mm的小骨料及其界面，是很少破裂的。基于大坝混凝土材料这些试验结果，如果将小骨料和固化水泥砂浆看作一种新的双相复合材料，即将小骨料和固化水泥砂浆等效为一种材料，并按分尺度分析方法预测多相复合材料的等效力学特性参数，使得大坝混凝土随机骨料模型中骨料的最小粒径增大到(20～40)mm，同时相应的骨料颗粒数目也大为减少，从而达到减少求解自由度的目的。应当指出，无论是“分尺度分析理论”或Eshel-by的“等效夹杂理论”，都只能预测弹性复合体的弹性模量和泊松比的等效弹性参数，不能给出其等效的强度。在细观力学的分析中，等效后新的复合基体的抗拉强度是关键参数之一