大坝混凝土动态性能细观力学及并行计算研究

1、114-Dec-21大坝混凝土动态性能细观大坝混凝土动态性能细观力学及并行计算研究力学及并行计算研究马怀发20062006年年6 6月月2424日日214-Dec-21一、研究的意义和试验成果分析二、混凝土细观力学研究方法三、混凝土试件数值模拟的基本方程四、混凝土试件的二维数值模拟五、混凝土三维细观有限元并行计算研究六、结束语汇汇 报报 大大 纲纲314-Dec-211 1、研究的目的和意义、研究的目的和意义1）重大工程抗震安全评价的主要环节 地震作用工程结构地震响应材料动态抗力 大坝混凝土动态抗力是大坝抗震安全评价的薄弱环节： 动态强度静态强度提高30%，依据国外早期少量湿筛试件的试验结果:

2、 大坝混凝土的湿筛试件和全级配试件动态性能存在差异。 2）拱坝受力特点： (1)高拱坝是一种空间结构，处于偏心受压状态。在强震时，坝体薄弱部位主要因坝体混凝土的动弯拉应力引起开裂，因此，动弯拉强度是高拱坝抗震设计中的控制性指标之一; (2)在地震作用时，坝体各部位都已经存在不同程度的静态应力。 因此,有必要研究: 在不同预静载作用下，混凝土动态性能的变化规律。414-Dec-211）小湾四级配混凝土和湿筛混凝土静动态弯拉试验成果：2 2、大坝混凝土的动态试验成果、大坝混凝土的动态试验成果514-Dec-212)小湾拱坝三级配混凝土，以冲击加载和三角波加载两种方式，在不同预静载水平下得到的动弯拉

3、强度：614-Dec-213）试验结果发现：a)全级配大坝混凝土的弯拉强度低于湿筛混凝土的弯拉强度,纯动弯拉强度增强系数都低于30%；b)全级配混凝土的动弯拉强度增强系数高于湿筛混凝土的动弯拉强度增强系数;c)一定预静载水平对动弯拉强度产生强化作用。4）数值摸拟的出发点：a)反映混凝土材料的内部结构特征，研究全级配混凝土与湿筛混凝土的细微结构差异；b)考虑混凝土材料的动载强化作用；c)反映预静载对混凝土试件动损伤破坏过程的影响。714-Dec-21一、研究的意义和试验成果分析二、混凝土细观力学研究方法三、混凝土试件数值模拟的基本方程四、混凝土试件的二维数值模拟五、混凝土三维细观有限元并行计算研

4、究六、结束语汇汇 报报 大大 纲纲814-Dec-21(1)微观层次， 材料的结构单元尺度在原子、分子量级，即从小于10-7 厘米至10-4厘米。(2)细观层次，尺度范围在10-4厘米至几厘米,或更大一些。在这个层次上，混凝土被认为是一种由粗骨料、硬化水泥砂浆和它们之间的过渡区（粘结带/界面）组成的三相材料。(3)宏观层次， 特征尺寸大于几厘米。材料假定为均质。要研究混凝土级配对混凝土宏观力学特性的影响，必须考虑混凝土细观结构的差异。水化硅酸钙混凝土的层次结构示意图微观宏观细观晶体原子结构混凝土颗粒 实验室尺寸 结构1 1、混凝土材料的细观结构层次、混凝土材料的细观结构层次914-Dec-21

5、2 2、细观数值模型、细观数值模型为了建立混凝土细微观结构特性与其在宏观力学特性的关系，自上个世纪70年代末，人们发展了混凝土细观力学研究方法。 细观力学将混凝土看作由粗骨料、硬化水泥胶体以及两者之间的界面粘结带组成的三相非均质复合材料，在细观层次上划分单元。考虑骨料单元、固化水泥砂浆单元及界面单元材料力学特性的不同，以其细观各相材料的破坏准则或损伤模型反映单元刚度的退化，利用数值方法模拟混凝土试件的裂缝扩展过程及破坏形态。1014-Dec-21(1)(1)二维随机骨料模型二维随机骨料模型瓦拉文用三维富勒级配曲线推导出了试件内截面任一点具有某种骨料直径的内截圆出现的概率，即瓦拉文公式。对于一个

6、特定的混凝土拌和物，使用这个公式即可产生试件截面上骨料的颗粒数。利用瓦拉文将三维问题简化为二维问题)0025.00045.0012.0053.0065.1()(10max1008max806max604max405.0max5.000DDDDDDDDDDPDDPkc骨料切面示意图AD0DxB1114-Dec-211)二维随机骨料模型多边形骨料模型 圆骨料模型1214-Dec-212）二维随机骨料模型细观有限元网格 骨料单元(灰色);砂浆单元(无色);界面单元(黑色)(b) 凸多边形骨料细观剖分(a) 圆骨料细观剖分1314-Dec-21(a)两级配(骨料含量47% )(b)三级配(骨料含量59

7、% )(c)四级配(骨料含量66%)随机球骨料分布模型(2)(2)三维随机骨料模型三维随机骨料模型1)随机球骨料模型分布将骨料简化为球体，按照实际级配的骨料粒径和各粒组骨料的比例，确定各种骨料的颗粒数。1414-Dec-212)随机凸多面体骨料模型分布在球形骨料模型的基础上形成多面体骨料基，以球骨料含量为控制参数，进在球形骨料模型的基础上形成多面体骨料基，以球骨料含量为控制参数，进行随机延凸，得到随机凸体骨料模型。行随机延凸，得到随机凸体骨料模型。(b) 三级配随机凸体骨料模型(a) 两级配随机凸体骨料模型随机凸体骨料模型1514-Dec-21三级配随机球骨料模型有限元网格图(a) 固化水泥砂

8、浆单元(b) 骨料单元(c) 界面单元(d) 对称切面网络图3)随机骨料细观有限元剖分1614-Dec-21(a) 骨料单元(b) 界面单元(c) 固化水泥砂浆单元(d) 试件两对称面上的细观单元剖面混凝土随机凸体骨料模型细观单元剖分1714-Dec-21一、研究的意义和试验成果分析二、混凝土细观力学研究方法三、混凝土试件数值模拟的基本方程四、混凝土试件的二维数值模拟五、混凝土三维细观有限元并行计算研究六、结束语汇汇 报报 大大 纲纲1814-Dec-21应变率（1/s）混凝土相对抗压强度与应变速率的关系相对抗压强度 相对抗拉强度应变率（1/s）混凝土相对抗拉强度与应变速率的关系1 1、混凝土

9、材料的应变率强化现象、混凝土材料的应变率强化现象（试验成果）1914-Dec-21 应变率效应是固体材料的基本特性。非均质材料较均质材料受应变率的影响更为显著。 混凝土在变形可逆的弹性阶段，几乎无应变率敏感现象。仅当加载至内部出现损伤时才呈现加载速率的影响。 应变率对混凝土动弹模的影响有类似的强化规律，但其影响弱于对动强度的影响。混凝土动拉压强度随应变速率增长，二者规律相似。2014-Dec-21混凝土抗拉强度和弹性模量与拉（或压）应变率的关系： (1)式中fts为混凝土材料的静拉强度，Es为静弹性模量。及表示混凝土动抗拉强度和动弹性模量是拉（或压）应变率的函数。Ht为强度强化系数； HE为弹

10、性模量强化系数，二者与应变率关系为： (2) 式中：At、Bt、Ct为强度强化参数，AE、BE、CE为弹性模量强化参数。2 2、混凝土的应变率强化定律、混凝土的应变率强化定律sEtsttEHEfHf)()()(lgexp)(lgexpEtCEEECtttBAHBAH2114-Dec-213 3、混凝土损伤演化定律、混凝土损伤演化定律应力 (MPa)轴向位移 / m混凝土轴拉试验（选自Schmidt-Hurtienne等人的文章）混凝土拉伸变形起控制作用2214-Dec-21采用了各向同性弹性损伤力学的本构关系来描述混凝土材料的力学性质。按照勒梅特应变等价原理，受损材料的本构关系： 式中E0为初

11、始弹性模量，为损伤后的弹性模量,D为损伤变量。 )1 (0DEEEuurrDmaxmaxmax0max0max00max11111100ftrrft应变双折线弹性损伤模型u应力（0D1）2314-Dec-214 4、混凝土试件数值模拟的动力学方程、混凝土试件数值模拟的动力学方程在动力学方程中，考虑到损伤弱化、应变率强化以及初始预静载对动态性能的影响。1）动力学微分形式： 其中)()()()()()(DPDPtPtUtKtUCtUMsdddddd )(tPd动载增量)( )()()(0ttUKDHDPdeeeeeeEd在动位移上因损伤产生的应力转移而产生的附加荷载 )()(0seeesUKDDP

12、在静位移上因损伤产生的应力转移而产生的附加荷载2414-Dec-212）动力学方程的积分形式： dVDdVDdSuTdVufdVuuuuDijVsklijkl)s(ijV)(klijkl(Si(t)iVi(t)ii(t)ii(t)iijVklijkl(t) 2514-Dec-21一、研究的意义和试验成果分析二、混凝土细观力学研究方法三、混凝土试件数值模拟的基本方程四、混凝土试件的二维数值模拟五、混凝土三维细观有限元并行计算研究六、结束语汇汇 报报 大大 纲纲2614-Dec-21计算模型（湿筛试件）弯拉试验及单元剖分示意图（单位:mm）5050150150150550150P1 1、预静载及加

13、载速率对混凝土动弯强度的影响、预静载及加载速率对混凝土动弯强度的影响2714-Dec-211 1）加载速率对动弯强度的影响）加载速率对动弯强度的影响动弯拉强度不但与加载速率有关，而且还与初始静荷载有关。 不同加载速率对动弹性模量和动强度的影响(a) 应力-位移曲线(b) 应力-位移局部放大曲线0.01.02.03.04.05.00.01.02.03.04.0位移(v/h 10-4)应力/ M Pa静载168kN/s200kN/s250kN/s300kN/s80%S+D250kN/s2.53.03.54.04.52.53.03.54.0位移(v/h 10-4)应力/ M Pa静载168kN/s2

14、00kN/s250kN/s300kN/s80%S+D250kN/s2814-Dec-212 2）初始预静载对动弯强度的影响）初始预静载对动弯强度的影响 计算表明：混凝土梁的动弯拉强度不仅与加载速率有关，而且与应变率、损伤积累密切相关。 预静载对动弯拉强度的影响（AE =0.12，At=0.20）(b) 应力位移局部放大曲线(a) 应力位移曲线0.01.02.03.04.05.00.01.02.03.04.0位移(v/h10-4)应力/ M Pa纯静载纯动载20%静载40%静载60%静载80%静载3.03.54.04.52.53.03.54.0位移(v/h10-4)应力/ M Pa纯静载纯动载2

15、0%静载40%静载60%静载80%静载100110120130140020406080100初始静载( % )相对弯拉强度( % )初始静预载对混凝土静动综合弯拉强度的影响2914-Dec-211 1）二维随机骨料随机参数模型）二维随机骨料随机参数模型(RAPM)(RAPM) (1)骨料随机分布（同随机骨料模型。按骨料级配及含量） (2)混凝土及其细观各相单元的抗拉强度和弹性模量随机分布这里所说的不均匀性包括混凝土骨料颗粒级配的差异和材料参数分布的统计随机性；随机骨料随机参数模型是在随机骨料模型的基础上，在细观有限元剖分后，再生成各类单元弹模和拉强度参数分布。 2 2、混凝土细观结构差异对其动

16、弯强度的影响、混凝土细观结构差异对其动弯强度的影响3014-Dec-212) 2) 细观结构对弯拉强度的影响细观结构对弯拉强度的影响 混凝土细观各相材料力学特性的不均匀性，使其动弯拉强度增强系数对初始静预载的敏感性增强。湿筛混凝土细观有限单元模型（单位: mm）细观不均匀性对静动综合弯拉强度的影响3114-Dec-213 3） 混凝土级配对动弯拉强度的影响混凝土级配对动弯拉强度的影响 全级配混凝土的动弯拉强度增强系数高于湿筛混凝土的动弯拉强度增强系数。 全级配混凝土骨料颗粒分布及单元剖分混凝土级配对静动综合弯拉强度的影响3214-Dec-214)4)骨料形态对混凝土对弯拉强度的影响骨料形态对混

17、凝土对弯拉强度的影响 1)湿筛混凝土试件 5050150150150550150P(a) 圆形骨料5050150150150550150P(b) 凸多边形骨料3314-Dec-21 湿筛混凝土两种细观骨料模型结算结果比较3414-Dec-212)全级配混凝土试件 四级配圆形随机骨料模型四级配凸多边形随机骨料模型3514-Dec-21 全级配混凝土两种细观骨料模型结算结果比较(b) 动弯拉强度增强系数随预静载变化曲线(a) 静动综合强度随预静载变化曲线3614-Dec-213 3、梁裂缝扩展分析、梁裂缝扩展分析1)1)静载作用下弯拉失稳过程静载作用下弯拉失稳过程(b) 加载到153.25kN(a

18、) 加载到152.25kN(c) 加载到153.75kN梁弯拉失稳破坏过程(d) 加载到极限荷载154.75kN3714-Dec-212)2)静动弯拉失稳破坏机理对比分析静动弯拉失稳破坏机理对比分析 (a)(a)静力静力: :作用时间相对较长，材料内部的损伤主要表作用时间相对较长，材料内部的损伤主要表现为微裂缝的扩展、联结和贯穿然后导致破坏；现为微裂缝的扩展、联结和贯穿然后导致破坏； (b)(b)动力动力: :因此材料内部的损伤主要表现为大量的微空因此材料内部的损伤主要表现为大量的微空洞和微裂纹生成，但来不及扩展。洞和微裂纹生成，但来不及扩展。 静载作用下梁的开裂过程动载作用时梁的开裂过程38

19、14-Dec-21一、研究的意义和试验成果分析二、混凝土细观力学研究方法三、混凝土试件数值模拟的基本方程四、混凝土试件的二维数值模拟五、混凝土三维细观有限元并行计算研究六、结束语汇汇 报报 大大 纲纲3914-Dec-211、并行计算研究的意义1）混凝土黏结带区域薄，一般在40-50微米；2）全级配混凝土试件尺寸大；全级配试件三维数值模型的节点自由度达到几十万甚至上百万，要进行大尺度全级配试件三维数值模拟，必须想办法解决计算速度这个瓶颈问题，必须进行并行计算研究。 4014-Dec-211、并行计算的两种方案1）对现有混凝土细观力学分析串行程序进行了并行算法优化设计（国防科大合作）：首先，从整

20、体上提出了一个将有限单元分布与未知量分布有机结合的首先，从整体上提出了一个将有限单元分布与未知量分布有机结合的整体并行算法设计方案整体并行算法设计方案。在进行并行算法设计时，整体上按有限元个数平均分配来进行任务划分，假设一共有p个处理器，将所有有限元分成p部分，每一部分中的有限元分配到一个处理器上。对一个给定的有限单元，如果分配到第k个处理器，则与之有关的计算在第k个处理器上进行。在对各个单元进行计算时，各个单元的计算之间是互不相关的，可以完全并行。然后分别针对刚度矩阵装配、双门槛不完全然后分别针对刚度矩阵装配、双门槛不完全Cholesky分解预条件、稀分解预条件、稀疏矩阵与向量相乘、稀疏向量

21、相加等核心算法，提出了相应的高疏矩阵与向量相乘、稀疏向量相加等核心算法，提出了相应的高效并行算法。效并行算法。模拟标准湿筛试件静加载损伤破坏过程，其三维模型共有模拟标准湿筛试件静加载损伤破坏过程，其三维模型共有71013节点，节点，78800单元，求解自由度约为单元，求解自由度约为21万。计算约需要万。计算约需要8个小时。计算平个小时。计算平台为水科究院网络中心台为水科究院网络中心Sun6800（4个主频个主频1G CPU、8G内存、内存、540G的磁盘整列存储空间）的磁盘整列存储空间） 4114-Dec-212）基于并行有限元程序自动生成系统 PFEPG （1）并行有限元程序自动生成系统 P

22、FEPG的基本思想： 元件化程序设计元件化程序设计 程序自动生成技术程序自动生成技术 区域分解法区域分解法将有限元程序分成：前处理分区、将有限元程序分成：前处理分区、start、bft、solv、E和和U等等 六六个模块个模块,其中其中E和和U元件程序根据用户给出的偏微分方程描述文元件程序根据用户给出的偏微分方程描述文件（件（vde或或pde），算法文件（），算法文件（nfe）由系统自动生成，其余几）由系统自动生成，其余几个元件程序由系统给出，不随表达式的变动而变动。个元件程序由系统给出，不随表达式的变动而变动。 4214-Dec-21（2）并行计算的程序结构 前处理数据主进程程序partit

23、ion元件bftm元件getpart元件msolv元件recv元件时间是否结束迭代是否结束从进程程序否sgetpart元件bfts元件U元件spart元件时间是否结束迭代是否结束否start元件E元件solv元件否否主调用程序图非线性动态问题并行求解流程框图4314-Dec-21（3）混凝土细观数值模拟的FEPG文件 gio命令流文件stdya 位移场对应单元计算文件为stdya.vdestdyb 应力场对应单元计算文件为stdyb.vde空一行#elemtype c8g2 单元类型，8节点六面体2阶高斯积分3dxyz 三维直角坐标系 4414-Dec-21 gcn命令流文件 defi 关键字

24、a stdya & a场求解位移，对应a场算法文件为stdya.nfe & 表示生成的文件不重新命名b stdyb b场求解应力，对应算b场法文件为stdyb.nfe空一行，表示定一段结束，以下通过命令流组织计算流程startilu a 对应预处理共轭梯度法不存单刚方式位移场初始化if exist stop del stop :1bft a 边值计算if exist end del end :2solvilu a 预处理共轭梯度法求解位移if not exist end goto 2 控制非线性迭代stress b 最小二乘法求解应力场post a 按时间步输出结果if not

25、 exist stop goto 1 控制时间循环4514-Dec-21 nfe算法文件 （1）位移算法stdya.nfe文件按照动静力学方程求解步骤，用FEPG有限元语言编写的求解位移的nfe文件。（2）应力算法stdyb.nfe文件在求得位移后，用最小二乘法求应力。 表示由已知变形算得的应力， 为待求的光滑后的应力。 dVudVVV)(0)(0u4614-Dec-21 vde文件 a)单元位移计算的stdya.vde文件b)单元应力计算的stdyb.vde文件根据上一时间步的应变率计算强度和弹模强化参数单元位移计算vde文件框图dfi信息段：定义未知量、坐标、函数以及传递已知量等对高斯点上

26、的应力、本构及损伤关系参量赋初值计算单元刚度矩阵和荷载列阵根据位移求高斯点主应力、主应变、判断该单元是否有新的损伤,更新最大应变计算损伤参数dfi信息段：定义未知量、坐标、函数、传递高斯点应力、最大拉应变及材料参数计算高斯点损伤参数根据高斯点上的本构及损伤关系参量计算单元刚度 单元应力计算vde文件框图根据高斯点应力和损伤参数计算最小二乘法的右端项将刚度阵赋值为零4714-Dec-215、PFEPG程序的生成和运行 1)填写Partition文件2)运行mpi命令：mpi gcn文件的前缀名3)把全部并行fortran源程序（*.f文件），fegen.b、runmpi、time0、nbefil

27、e、partition.dat、*.io传送到并行机4)在并行机上运行命令: csh x fegen.b5)运行前处理6)修改hosts文件和runmpi文件7)在并行机上运行如下命令: csh x runmpi8)运行后处理4814-Dec-21(1)并行计算环境 1）硬件环境： 4台以主频3G、内存1G 的PC机以千兆以太网连接的工作站。 2）软件环境： Linux操作系统,PFEPG并行有限元程序自动生成系统平台采用MPICH通信。6、并行计算分析4914-Dec-21(2)细观模型选取及并行程序生成 1）计算模型：整个剖分区域共有111988个单元，76810个节点。 图55 三维随机

28、骨料分布湿筛混凝土弯拉三维数值模型(c) 空间网格(d) 对称面网格切面(f) 骨料单元(e) 固化水泥砂浆单元(g) 界面单元P(a) 三维数值模型平面图（单位: mm）5050150150150550150(b) 约束条件5014-Dec-21a）一种单元类型对应的一个vde文件；b）跨中细观剖分区域与其它区域独立分开；c）执行gio命令流文件后修改pre文件和对应各物理场以io为扩展名的文件；d）再执行pre命令和其它vde文件（除gio文件中vde文件）e）生成单机串行程序，再按PFEPG程序步骤生成该问题的并行程序。f）修改材料参数（修改elem0和elemb0中的材料参数），进行并行计算。 5114-Dec-21模型并行计算分区图partition.dat文件为4 2 4个分区，2表示按坐标分区2 2 1 x方向剖分数为2，y方向剖分数为2，z方向剖分数为15214-Dec-21(3)预静载作用下三维细观模型并行计算分析 三维模型并行计算是在主频为3G、内存为1G 的4个PC节点的并行集群上进行的。最大控制位移为5.010-4，非线性迭代位移误差范数控制为10-15，最大迭代次数为15次。静力计算耗时约15.6小时；动力计