混凝土梁四点弯曲试验的并发多尺度区域分解法模拟

邱莉婷，马福恒，沈振中，张 湛，霍吉祥

(1.南京水利科学研究院，江苏 南京 210029； 2.河海大学 水利水电学院，江苏 南京 210098； 3.河南省西霞院水利枢纽输水及灌区工程建设管理局，河南 郑州 450008)

混凝土作为由硬化水泥砂浆、骨料及二者界面过渡区组成的非均质三相复合材料，在浇筑过程中由于温度应力和收缩问题会形成大量的孔隙和原始微裂纹[1]。其断裂是在外界荷载作用下，局部微缺陷发展成核、扩展汇聚，宏观裂纹形成，裂纹失稳扩展的交织发展过程。整个断裂过程无大比例塑性流动产生，仅吸收了裂缝表面形成过程中的剩余能量，是典型的准脆性断裂[2]。准脆性混凝土材料由于细观结构的多相性和不均匀性，其宏观力学响应十分复杂，表现为非线性、应变软化和应变局部化等材料特性。其非线性力学特性主要有弹塑性、黏性、损伤和断裂等几种，损伤占有重要地位。损伤力学是研究材料的微、细观缺陷及裂纹萌生和发展对其宏观力学特性的影响，以及材料和结构损伤演化过程及机理的破坏力学。当有限元计算采用局部损伤模型时会存在网格敏感性和零能耗问题[3]。隐式梯度损伤模型[4]属于非局部损伤模型，以变形的梯度形式反映微结构变化和相互作用，把空间相互作用与梯度公式的计算效率相结合，适用于同非线性有限元法结合开展混凝土损伤失效的数值模拟。

同时，传统单一网格的有限元模型无法对混凝土线弹性区域和损伤集中区域进行区别处理，不能将有限的计算资源用于需要重点分析的局部易损区域。而采用过渡单元进行线弹性区域粗网格和非线性区域细网格的连接，当网格尺寸跨度较大时，过渡单元的网格剖分质量难以保证。如进行自适应有限元分析[5]往往需要进行多次自适应网格剖分才能达到要求的计算精度。对此，并发多尺度方法[6]可以使高效率的粗网格和高精度的细网格在同一有限元模型中共存，且对混凝土的应变软化行为不存在层级多尺度方法[7]所面临的尺寸依赖性问题，可用于混凝土的损伤失效全过程分析。但计算过程中直接包含细网格子区域仍会使计算资源消耗较大。整体有限元撕裂对接法（Total Finite element tearing and interconnecting method,Total-FETI）[8]属于局部边界条件不重合的区域分解法，基于分而治之思想，具有高度的并行性和并行效率，适合于大规模数值计算。

本文采用隐式梯度损伤模型描述混凝土材料的非线性本构关系，将并发多尺度方法和整体有限元撕裂对接法结合，在兼顾计算精度和计算效率的基础上，构建混凝土损伤失效分析的并发多尺度区域分解模型，并将该模型应用于混凝土预制切口四点弯曲梁试验的数值模拟。

1 混凝土单边切口梁四点弯曲试验数值模型

1.1 隐式梯度损伤模型在Total-FETI中的实现

有限元撕裂对接法[9]将所研究区域分解为多个不重叠的独立子区域，通过拉格朗日乘子去除各子区域界面间的连续性约束。其求解由两部分组成：首先是子区域界面问题求解，即采用预条件共轭梯度算法（Preconditioned conjugate gradiant，PCG）求解仅含界面未知量的并由狄利克雷边界条件约束的缩减系统；然后再进行由纽曼边界条件约束的各个独立子区域问题求解。对于界面问题，在PCG迭代求解过程将误差进行全局化处理，从而加快计算收敛速度。同时，相比传统的子结构法，FETI的子区域划分不需要人为干预，也不要求满足子结构内部节点先编号、界面节点后编号的编号原则，无需再进行节点编号排序优化，且子区域间通信规模也较小。FETI法的内在机制使得计算时的收敛快慢和子域划分的数量多少无直接联系，可以根据计算需求灵活地进行子区域划定。进一步地，整体有限元撕裂对接法（Total-FETI）[8]在有限元撕裂对接法基础上，将拉格朗日算子不仅用于子区域交界面的约束，同时也用于边界子区域的狄利克雷边界条件施加，从而达到简化子区域刚度矩阵的求逆过程，这也是本文多尺度区域分解算法构建的基础。

本文主要研究混凝土损伤失效的非线性演变过程，因此，这里给出非线性Total-FETI离散后的牛顿拉弗森迭代（Newton-raphson）求解公式。假定在某加载步的第k+1牛顿拉弗森迭代步，其位移向量和拉格朗日乘子的更新如下：

式中：u为位移向量；拉格朗日算子λ同时考虑了子区域界面位移自由度和狄利克雷边界条件的施加。

非线性有限元撕裂对接法的线性方程系统和子区域Ω(s)间的位移连续条件可表述为：

式中：A(s)为切线刚度矩阵；δu(s)和δλ分别为位移增量和拉格朗日乘子增量；布尔矩阵B(s)由子区域界面的布尔矩阵和狄利克雷边界条件对应的布尔矩阵按行串连组成；fext(s)为外力向量；fint(s)为内力向量；Ns为子区域数量。

隐式梯度损伤模型属于特殊的非局部损伤模型，其求解非局部等效应变的隐式梯度方程将非局部模型中的权函数替换成格林（Green）函数，也称修正的Helmholtz’s方程[10]。

式中：非局部等效应变ɛeq不是由等效局部应变ɛ及其导数显式描述，而是作为包含方程（4）的边值问题和近似边界条件的解；梯度参数为内部长度参数l的函数，其取值为c(l)=l2/2。

Total-FETI由于考虑了子区域界面相容条件，其控制方程的变分形式以混合变分形式表述。Total-FETI与隐式梯度损伤模型结合后，由于子区域界面的非局部等效应变ɛeq的存在，需要对给定子区域界面的拉格朗日算子进行如下修正：

式中： λεeq考虑了非局部等效应变的自由度。

通过扩展布尔矩阵B(s)将λd装配到合适位置，布尔矩阵Bd(s)由原来Total-FETI的子区域界面的布尔矩阵B(s)和非局部等效应变对应的布尔矩阵Beq(s)按行串连组成：

1.2 并发多尺度方法与Total-FETI的结合

细网格子区域和粗网格子区域在其连接界面需满足位移连续条件。细网格子区域界面节点位移可以通过粗网格子区域界面单元的形函数和界面节点位移进行插值求解。这里将粗网格子区域界面节点称为主节点，将细网格子区域界面未能与主节点匹配的其余节点称为从节点。从节点的自由度由线性多点约束以不同的插值方式实现与主节点的耦合。其中主节点间的约束为线性多点约束，从节点和主节点之间的约束称为尺度间的线性多点约束。对于两个有限元网格尺寸不一致的子区域界面，其尺度间的线性多点约束用矩阵形式[11]可表示如下：

式中：u(s)为子区域 Ω(s)的位移向量；矩阵P(s)由线性多点约束组成，建立了从节点和主节点间的联系，即细网格节点自由度通过粗网格子区域交界面处单元的形函数插值求得。

尺度间的线性多点约束可通过拉格朗日乘子法实现[12]。采用拉格朗日乘子在Total-FETI中添加额外数量方程，通过定义一个修正的布尔矩阵来实现，修正后布尔矩阵由Bd(s)和约束矩阵P(s)按行串连组成：

扩展后的拉格朗日乘子U包含邻域间主节点的界面约束λd及主节点与从节点间的约束η：

在引入隐式梯度损伤模型后，Total-FETI的界面问题求解需要重新修正。隐式梯度损伤模型在每个节点需要额外的自由度来描述其非局部等效应变，所以子区域的刚度矩阵是非对称矩阵。从而柔度矩阵和子区域界面问题的拉格朗日矩阵均为非对称矩阵，对其进行迭代求解时只有双共轭梯度稳定解法（BiCGStab）和广义极小剩余法（GMRES）两种迭代求解器可满足计算要求。考虑到子区域界面方程的迭代求解难以构建合适预处理算子，在整体有限元撕裂对接法求解框架下，不妨回到有限元撕裂对接法的最初形式：

式（12）可以通过串联不同子区域的矩阵和向量改写成熟悉的有限元支配方程形式Au=f。方程组采用分块矩阵的形式表述如下：

将各个子区域的刚度方程和描述各子区域界面和边界条件的布尔矩阵直接组装成有限元支配方程的整体系数矩阵，称为双重组装法[13]。再采用非对称稀疏线性方程组的并行直接求解器Pardiso[14]进行大型线性方程组的求解。采用双重组装法的整体有限元撕裂对接法不再需要进行子区域的界面问题和刚体位移求解。

1.3 并发多尺度区域分解算法建立流程

混凝土损伤分析的并发多尺度算法建立流程主要包括子区域分解及多尺度有限元网格剖分，子区域界面信息和狄利克雷边界信息提取，布尔矩阵及其转置矩阵建立及各加载步的牛顿拉弗森迭代求解。并发多尺度有限元网格剖分，首先按照构件几何特征和重点分析区域进行子区域分解，其子区域的几何形状可以是不规则的；再对全部子区域进行粗网格剖分，粗网格的剖分要求能较好地反映构件几何及子区域边界形状，也要能准确捕捉构件线弹性阶段位移场的空间分布和变形梯度变化；重点分析区域在保留粗网格子区域边界节点的基础上进行指定精度的细网格剖分，细网格子区域的内部可以按照计算精度要求进行不同类型的有限元网格剖分，不受原粗网格子区域内部节点限制。

2 计算分析

为验证模型的合理性，在前人研究[15-16]的基础上，对图1所示的混凝土单边切口梁四点弯曲试验进行数值模拟验证。混凝土试件尺寸为500 mm×100 mm（长×高），净跨为450 mm，底部预制切口高度d为10 mm。模型中预制切口采用去除节点水平向约束的零厚度切口模拟。试件底部右侧支撑处和左侧支撑处节点均施加竖直向和水平向约束。数值计算过程采用位移加载方式控制，在梁顶部两个三分点处的节点分别施加强制位移荷载。为与文献[16]作对比分析，采用试件底部预制切口偏左侧7.5 mm处节点的竖向位移作为梁的挠度。

图1 混凝土单边切口梁四点弯曲试验的几何及边界条件示意（单位：mm）Fig.1 Geometry and boundary conditions for four-point bending test (unit: mm)

首先，进行并发多尺度区域分解法建模。如图2(a)所示，将试件分解为10个独立的子区域，其中试件中部的子区域数量较多，主要是有限元网格精度区分需要，以及验证子区域分解是否影响模型计算结果。各子区域根据计算精度要求进行不同尺寸的有限元网格剖分，剖分情况如下：左、右两侧的子区域Ω(1)、Ω(10)均采用h=10 mm的有限元网格；子区域 Ω(2)、Ω(5)、Ω(6)、Ω(9)均采用h=5 mm 的有限元网格；梁中部可能发生应变局部化的4个方形子区域Ω(3)、Ω(4)、Ω(7)、Ω(8)的网格尺寸一致，为验证模型是否存在网格敏感性问题，将这4个子区域分别剖分为5.000、1.250 和0.625 mm共3种不同精度的有限元网格进行对比分析。支撑点和加载点均采用宽为20 mm，顶部高出梁10 mm的刚体进行模拟，其弹性模量取值比混凝土梁的要大一个数量级，支撑点处的网格进行适当局部加密处理。本算例假设为平面应力问题进行计算，各子区域有限元计算网格均采用四节点的四边形单元。为直观地说明有限元网格剖分情况，在图2(b)中给出了梁中部子区域网格为1.250 mm时的有限元网格剖分示意。

图2 混凝土单边切口梁子区域分解及有限元网格剖分示意（单位：mm）Fig.2 Domain decompostion and finite element of four-point bending beam (unit: mm)

隐式梯度损伤模型计算所采用的参数说明[16]如下：采用修正的Von Mises模型计算等效应变，其混凝土的压缩和拉伸强度比值η为10；采用Mazars模型中的损伤演化法则，其控制软化速率的β取值为300，控制残余应力的α取值为0.92，损伤发生时的应变阈值k0取值为0.000 075；弹性模量E取值为40 000 MPa，泊松比ʋ为0.2，梯度参数c=4 mm2。双重组装后线性方程组系统的求解采用直接求解法进行求解，求解器为Pardiso。

计算得知，各加载步的牛顿拉弗森迭代步数均在4个加载步以内计算收敛，可知模型鲁棒性好。梁中部3种不同精度有限元网格剖分模型所得的损伤分布差别不明显，这里仅将中等精度网格h=1.250 mm的损伤分布和发展过程在图3中展示。由图3可知，损伤萌生于预制切口的尖端，并迅速沿梁中轴线向梁顶部扩展。由于试件结构和加载方式的对称性，损伤场关于梁中轴线对称分布，损伤模式与Simone等[16]的模拟结果（见图4，图例为损伤变量，0为未损伤，1为完全损伤）基本一致。同时，损伤分布在子区域界面连续性良好，未表现出子区域分解敏感性。

图3 混凝土单边切口梁四点弯曲试验的损伤扩展过程（h=1.250 mm）Fig.3 Damage evolution of four-point bending beam (h=1.250 mm)

图4 Simone文献混凝土单边切口梁四点弯曲试验的损伤扩展过程（h=5.000 mm）[16]Fig.4 Damage evolution of four-point bending test in the paper of Simone (h=5.000 mm)[16]

图5所示非局部等效应变分布和发展过程表明，应变集中区域与损伤集中区域一致，均位于梁纯弯段的底部预制切口尖端区域。且非局部应变的分布在子区域界面保持完好连续性，受子区域分解的影响较小。

图5 混凝土单边切口梁四点弯曲试验的非局部等效应变发展（h=1.250 mm）Fig.5 Nolacal equivalent strain evolution of four-point bending beam (h=1.250 mm)

梁中部选用不同精度网格模型所得到的荷载-竖向位移曲线分别如图6中的实线所示。峰值荷载前的线性增长阶段，试验结果与数值模拟结果的一致性较好；对于峰值荷载，其中试验峰值荷载为3.68 kN；文献数值模拟峰值荷载为3.96 kN；本文h=5.000 mm峰值荷载为3.81 kN，h=1.250 mm峰值荷载为3.60 kN，h=0.625 mm峰值荷载为3.50 kN。可知峰值荷载随着网格精细程度的增加有小幅减小，但与试验相差2.2%～4.9%；软化阶段h=5.000 mm与参考文献的数值模拟结果相近，h=1.250 mm和h=0.625 mm的软化曲线几乎重合，试验软化曲线介于h=5.000 mm和h=0.625 mm的软化曲线之间。综上可知，基于隐式梯度损伤模型的并发多尺度区域分解法不存在网格敏感性，且计算结果受子区域分解影响有限，这表明模型是合理的。

图6 混凝土单边切口梁四点弯曲试验的荷载-竖向位移曲线Fig.6 Load-deflection curves of four-point bending beam

3 结 语

（1）采用隐式梯度损伤模型描述混凝土材料的非线性本构关系，并在易损伤的局部子区域预设高精度有限元网格，同时，采用尺度间线性多点约束进行不同网格精度子区域的界面协调连接，构建了混凝土损伤分析的并发多尺度区域分解模型。

（2）将模型用于混凝土单边切口梁的四点弯曲试验模拟，对可能损伤的局部子区域分别采用了3种不同精度的有限元网格计算，结果表明模型不具网格敏感性。同时，损伤及非局部等效应变分布受子区域分解影响不显著，在子区域界面的连续性良好，荷载-竖向位移曲线的计算结果与数值模拟结果的一致性较好。计算结果表明模型是合理的，可用于混凝土损伤分析的并发多尺度区域分解计算。

（3）本文建立的并发多尺度区域分解法基于局部易损子区域的高精度有限元网格预设。考虑到混凝土材料细观结构的高度非均质性，其裂缝扩展路径具有不确定性，数值计算过程中需要不断地进行网格自适应更新，以满足损伤局部化区域的高精度网格剖分要求。因此，建立子区域的自适应更新方法，实现损伤子区域高精度有限元网格的实时自适应更新是下一步研究工作的重点。