0 引 言

混凝土作为水利工程的常用材料，在温度、地震等作用下常常产生裂缝。大部分的裂缝并不影响建筑物的正常使用，而少部分裂缝可能会影响结构的使用寿命。研究混凝土的断裂机理，是目前比较热门的课题。现有的混凝土断裂力学认为：混凝土的断裂特性表现为具有时效历程的断裂过程[1]。针对混凝土具有的准脆性的特性，Hillerborg[2]建立了描述混凝土断裂过程区传递应力与裂缝宽度关系的虚拟裂缝模型，该模型以一条软化曲线来描述混凝土断裂过程区的软化特性，成为后来学者研究混凝土断裂的常用工具[3-5]。国内外确定软化曲线的方法主要有：直接拉伸试验法、J 积分法及逆推分析法。其中 Evans[6]等最早采用直接拉伸试验获得了混凝土拉伸应力-应变曲线的下降段，林皋[7,8]等研究了不同应变速率下直拉试验中混凝土的力学响应，Rice[9]等在考虑裂缝端塑性变形的情况下应用全量理论提出了求解平面裂缝问题软化曲线的 J 积分方法。其后的学者研究发现，J 积分方法获得的软化曲线准确性较差，而采用直接拉伸法时由于混凝土的非均质特性性，容易面临应力集中、试验机刚度不足等问题导致难以获得软化曲线完整的下降段[10]，因此 Roelfstra[11] 提出通过反复调整软化曲线参数以使数值仿真得到的荷载位移曲线同试验测得的荷载位移曲线之间差异最小的方法来确定混凝土的软化曲线，称为逆推法。国内学者也对此进行了较为有益的探索[12,13]。由于混凝土断裂问题本身具有的复杂性，目前的研究较难取得统一的认识。因而，深入研究混凝土的断裂机理有其必要性。

 

本文以双线性软化曲线为例，通过 Python 脚本自编程，直接在 Abaqus 平台上实现了双线性软化曲线控制参数的逆推寻优计算。并将逆推得到的荷载-裂缝口张开位移曲线 P-CMOD 及断裂能同试验测得的数据进行了比较分析。

 

1 双线性软化曲线及控制参数

与传统的线弹性断裂力学理论不同，非线性断裂力学理论认为混凝土开裂破坏时裂缝尖端存在具有一定几何尺寸的微裂区，称之为断裂过程区（FPZ）。断裂过程区介于完全丧失黏聚力的宏观裂缝区与未开裂区混凝土之间，其力学特性表现为随着微裂区持续张开断裂过程区的黏聚应力持续降低。混凝土断裂过程区的存在使得混凝土实际断裂过程中存在明显的软化特性。作为描述混凝土断裂过程裂缝尖端软化本构关系的软化曲线，常有几种简化形式主要包含直线型、折线型、指数型软化曲线。其中折线型软化曲线中的双线性软化曲线是混凝土断裂分析广泛使用的本构模型，如图 1 所示，其数学表达式为：

 

 

 

这种形式相应的断裂能可以表示为：

 

 

 

可以看到双线性软化曲线的纵、横坐标截距、及转折点处的位置坐标四个参数是确定软化曲线形状及断裂能大小的基本参数。

 

 

图 1 双线性软化曲线

 

2.1 误差定义

双线性软化曲线参数不同，代入有限元程序进行数值计算获得的 P-CMOD 曲线就不同。逆推分析的目的即为找出双线性软化曲线最合适参数，以使数值计算得到的 P-CMOD曲线能够同试验测得的 P-CMOD 曲线之间误差最小。为了刻画它们两者之间的差异，记数

 

 

图 2 误差定义

 

2.2 优化算法

软化曲线的逆推问题可简化为找寻一组参数 ，以使得 最小的单目标优化问题，其数学表达式为：

 

 

 

 

式中， 为常数。在优化过程中， 越小，个体的适应度越大，适应度越大的个体在选择过程中越容易得到保留。

 

2.3 逆推流程

本文采用 Python 编程，将所编脚本直接在 Abaqus 平台上运行计算。构建的逆推运行流程如下：①根据三线性软化曲线控制参数确定优化变量；②配置进化算法相关参数（杂交、变异、选择算子、终止条件等），生成初代种群并对个体进行编码；③根据该代种群相应个体确定软化曲线代入有限元程序中进行计算得到不同个体数值仿真 P-CMOD 曲线与试验 P-CMOD 曲线的误差 ；④根据种群中个体误差 计算该代的适应度；⑤根据适应度分布状况进行杂交、变异运算、选择运算生成下一代种群并计算相应适应度，⑥判断该代优化变量是否满足终止条件（达到最大进化代数或 达到允许误差）：若不满足，则转到流程③；若满足，对最后一代种群中的最优个体进行解码，即为软化曲线的最优参数。

 

3 数值算例

3.1 模型建立

逆推的实质是比较数值计算得到的 P-CMOD 曲线同试验测得的 P-CMOD 曲线的差

异，本文选取 Jeffrey[15]的梁深为 63mm 的三点弯曲切口梁试验作为数值算例，编号为

CB63-80a、CB63-80b、CB63-80c 的三个试件梁长 L 为 350mm，跨长 S 为 250mm，梁深D 为 63mm，梁厚 t 为 80mm，跨中切口深α0为 21mm，试件尺寸参数如图 3 所示。试验测得的 P-CMOD 曲线见图 4，根据 RILEM 推荐的方法计算的断裂能分别为 123 N/m、124 N/m、123 N/m。选取黏聚裂缝模型作为逆推分析中的计算模型，其 值根据专著[16]介绍依据 Brazilian test 确定，取 。计算过程中混凝土其余部分取弹模，泊松比 。逆推分析过程中定义误差取 ，构造适应度函数取 。

 

 

图 3 三点弯曲梁试件尺寸

 

 

图 4 梁深 D=63mm 试件的 P-CMOD 试验曲线

 

3.2 逆推结果及讨论

基于 2.3 节所阐述的流程进行逆推分析，其结果如表 1 所示，CB63-80a、CB63-80b、CB63-80c 三个试件计算的适应度分别为 0.973、0.975、0.934。图 5、图 6、图 7 分别为逆推得到的 P-CMOD 曲线同试验测得的 P-CMOD 曲线的比较，可以看到，两者均能够较好吻合，试验测得的荷载峰值分别为 2717.76N、2730.73N、2822.60N，逆推得到的荷载峰值分别为 2644.14N、2737.61N、2683.32N，两者之间的误差均分别为 2.7%、0.25%、4.9%。逆推分析的断裂能可以根据表 1 得到的软化曲线控制参数利用式（2）计算得到，分别为 134.7N/m、124.3N/m、124.9N/m，文献[15]中根据 RILEM 推荐方法计算的断裂能分别为 123 N/m、124 N/m、123 N/m，两者之间的误差分别为 9.5%、0.24%、1.5%，见图8。逆推得到的 P-CMOD 曲线的荷载峰值同试验测得的数据之间误差均不超过 5%，逆推得到的断裂能同试验计算的断裂能之间误差不超过 10%。Kyoungsoo[17]根据 Jeffery[15]的试验数据研究了软化曲线转折点的确定，其首先根据 Shah[18,19]的两参数断裂模型(TPFM)计算了断裂能，接着依据 Wittmann

[20]的研究确定软化曲线转折点处 为 0.25，进而确定软化曲线的形状代入有限元程序中进行了数值计算，最后比较了该方法数值计算得到的P-CMOD 无量纲曲线（将原 P-CMOD 曲线原横纵坐标无量纲化处理，处理后的横坐标表示为 ，纵坐标表示为 ）同试验测得的曲线之间的差异，图 9 比较了采用逆推法得到的 P-CMOD 无量纲曲线同 Kyoungsoo 文献[17]的方法得到的 P-CMOD 无量纲曲线之间的差异，可以发现逆推得到的 P-CMOD 无量纲曲线比文献[17]采用的方法得到的 P-CMOD 无量纲曲线更接近试验测得的实际曲线，体现了逆推法的优势。图 10 比较了逆推得到的三条软化曲线，以转折点处为区分可将双线性软化曲线划分为第一下降段和第二下降段，图 10 的结果显示，逆推得到的三条软化曲线的其第一下降段相比十分接近，Ba?ant 的专著[16]中认为双线性软化曲线的第一下降段的斜率决定了数值计算的 P- CMOD 曲线的荷载峰值，本文的计算结果也从侧面证了这点。

 

表 1 逆推分析结果

 

 

 

图 5 CB63-80a 试件逆推 P-CMOD 曲线同试验曲线比较

 

 

图 6 CB63-80b 试件逆推 P-CMOD 曲线同试验曲线比较

 

 

图 7 CB63-80c 试件逆推 P-CMOD 曲线同试验曲线比较

 

 

图 8 逆推断裂能与 RILEM 方法计算的断裂能之间的比较

 

 

图 9 试验与逆推得到的 P-CMOD 无量纲

 

 

图 10 CB63-80a、CB63-80b、CB63-80c 试件

逆推软化曲线之间的比较

 

4 结论

本文采用 Python 编程利用进化算法对混凝土双线性软化曲线进行了逆推分析，可得出如下结论：

（1）与常规的直接编程方法相比，本文采用 Python 脚本编程，直接在 Abaqus 平台上实现逆推分析计算，避开了繁琐的有限元求解器的编程工作，专注于参数化建模的实现，提高了分析计算的效率，是一种值得推广的方法。

（2）计算结果表明，采用本文的方法对混凝土双线性软化曲线进行逆推分析是行之有效的，逆推得到的 P-CMOD 曲线荷载峰值、断裂能同试验测得的荷载峰值、断裂能之间误差不超过 10%，满足工程应用的实际要求。

（3）利用 Abaqus 丰富的二次开发接口，能够实现例如定义本构、编制算法等功能。继续深入研究混凝土的断裂本构关系、更精细地定义误差以及改进优化算法等是下一步研究的重点。

 

资料来源：达索官方