基于串联均质化模型的再生混凝土抗压强度和破坏机理分析

Analysis of Compressive Strength and Failure Mechanism Based on Series Homogenization Model for Recycled Concrete

摘要: 本文针对再生混凝土材料，利用均质化分析方法，建立了串联均质化模型，采用新型有限元法——基于余能原理的基面力单元法对再生混凝土进行细观层次的损伤研究，利用所编非线性损伤基面力单元法分析程序，分析了再生混凝土在单轴压缩荷载情况下的应力应变软化曲线和破坏机理，并将计算结果与试验和细观随机骨料模型分析方法的结果进行对比分析。研究结果表明：这种串联均质化分析方法可以用于分析再生混凝土材料的非线性损伤分析，可以得到再生混凝土材料的抗压强度、应力应变软化曲线和裂纹扩展的过程。

Abstract: The homogenization analysis method was used to establish a series homogenization model for recycled concrete materials in the paper. A new type of finite element method—Base Force Element Method based on complementary energy principle is used to conduct meso-level damage research on recycled concrete. The stress-strain softening curve and failure mechanism of recycled concrete under uniaxial compression load are analyzed using the nonlinear damage analysis program of base force element method. The calculation results are compared with the results of experiments and meso-level random aggregate model analysis methods. The research results show that this series homogenization analysis method can be used to analyze the nonlinear damage analysis of recycled concrete materials. The compressive strength, stress-strain softening curve and crack propagation process of recycled concrete materials can be obtained using the present method.

1. 引言

再生骨料混凝土(Recycled Aggregate Concrete, RAC)简称再生混凝土(Recycled Concrete)，是将废弃混凝土经过清洗、破碎、分级和按一定比例与级配混合形成再生骨料，部分或者全部代替砂石等天然骨料配制成的新混凝土。它作为一种绿色环保型建筑材料已经得到广泛的重视 [1]。

在再生混凝土性能研究方面，国内外学者对其展开了大量的试验研究，并取得了一些研究成果 [1]。但试验周期长，费用高，材料内部应力、应变和破坏机理测量较难。因此，开展再生混凝土数值模拟研究是十分有益的工作。2009年，肖建庄 [2] 建立了再生混凝土随机集料模型，采用有限元法模拟得到了再生混凝土在单轴受压试验下的应力–应变曲线。2012年，李文贵 [3] 通过变参数数值模拟分析研究了再生混凝土五相介质的力学参数对模型的应力分布特征的影响。彭一江等 [4] [5] 推导出基面力元法并采用这种新型有限元法对再生混凝土进行了数值模拟分析，分别在静力与动力荷载作用下进行单轴拉压试验，研究其力学性能。2015年，Liu [6] 在细观尺度下，研究了采用不同形状骨料再生混凝土的受压性能。2016年，Rajendra [7] 建立了虚拟裂缝模型和双K断裂模型，确定了再生混凝土不同粗骨料含量的断裂参数。2018年，谢永滨 [8] 对再生混凝土棱柱体进行了大量试验分析不同取代率下的再生混凝土应力–应变关系。Anuruddha [9] 通过数值模拟分析得到再生混凝土的弹性模量和强度随老砂浆的含量增加而降低。Job et al. [10] 通过数值模拟及回归分析的方法，研究了再生混凝土的力学性能与再生骨料取代率的关系。2019年，Tan [11] 基于离散元方法，对再生混凝土进行了二维数值模拟，主要研究了新老界面薄弱环节对再生混凝土损伤破坏过程的影响。Guo [12] 针对再生混凝土的非均质性建立了徐变耦合模型，研究再生骨料对再生混凝土的徐变的影响。F. Kazemian et al. [13] 对再生混凝土进行了抗压强度、抗弯强度和断裂能试验研究。对比了经过处理和未经过处理的再生混凝土的力学性能。2021年，Alireza [14] 采用离散元法建立了再生混凝土数值模型，使用研究了该模型在无侧限的压缩试验下和加卸载条件下的宏观与微观力学特性。还有一些学者在这一领域开展研究工作。

本文拟针对再生混凝土材料，利用均质化分析方法，从细观层次建立串联均质化模型，运用一种基于余能原理的基面力单元法 [4] 进行再生混凝土抗压强度和破坏机理研究。

2. 再生混凝土串联均质化模型

2.1. 串联均质化模型

基于细观损伤单元等效化模型的基本思想，本文采用简单的材料强度等效方法求解Reuss串联模型，将再生混凝土非均质的单元进行均质化，建立细观损伤等效模型，如图1所示，其中E为弹性模量， $\sigma ,\epsilon$ 分别为应力与应变。

Figure 1. Reuss series model

图1. Reuss串联模型

由本构方程 $\epsilon =\frac{\sigma }{E}$，基于等应力假设，忽略横向变形可以得到

$\epsilon =c_1{\epsilon }_1+c_2{\epsilon }_2$ (1)

$\frac{\sigma }{E}=c_1\frac{\sigma }{E_1}+c_2\frac{\sigma }{E_2}$ (2)

$\frac{1}{E}=\frac{c_1}{E_1}+\frac{c_2}{E_2}$ (其中 $c_1=\frac{V_1}{V}$， $c_2=\frac{V_2}{V}$ ) (3)

具体实施时，先进行大网格划分，再进行细小网格划分，然后将大网格进行均质化处理，如图2所示。粗网格的剖分方法与判断规则与细网格的相同，判断的规则是根据单元节点位置确定的，当单元四个节点中有大于等于3个节点投影在骨料(或老砂浆，或新砂浆)介质时，那么判定该单元属性为骨料(或老砂浆，或新砂浆)；当单元节点有落在骨料介质中也有落在老砂浆中时，该单元定义为老界面单元；同理当单元节点有落在老砂浆介质，又有节点落在新砂浆介质时，该单元便定义为新界面。判定细网格单元属性之后，统计各属性的细小网格数量，根据网格尺寸计算面积量，并计算各组分在粗大网格的占据比例。采用均质化方法将粗大网格等效成均匀的单一属性单元。

采用FORTRAN语言编程，将划分网格与单元属性的赋予编写成程序进行计算，得到试件全部的网格节点编号与坐标，以及单元属性文件由计算程序可以得到单元的属性分布，并得到含有多相介质非均匀单元的组成，为后续的均质化计算提供数据。

Figure 2. Mesoscopic model of regenerated concrete meshes

图2. 再生混凝土网格划分细观模型

2.2. 损伤本构模型

本文对再生混凝土进行数值模拟，由于再生混凝土材料的非均匀性，在单轴受压时，应力趋近于峰值时表现出极强的非线性，因此本文采用多折线应力应变本构关系进行计算。引入标量损伤变量D，由于拉伸或者压缩导致材料受损，受损材料的柯西应力与有效应变的关系为

$\sigma =E_0\left(1-D\right)\epsilon$ (4)

若忽略损伤对泊松比的影响，损伤后的弹性模量可以用初始弹性模量表示

$E=E_0\left(1-D\right)$ (5)

式中， $E_0$ 表示初始弹模，E为损伤后弹模。因此对于再生混凝土中的五相介质的损伤弹性模量可以表示为(新砂浆(m)、骨料(ag)、老界面(oitz)、老砂浆(om)、新界面(itz))。

$\{\begin{array}{l}{E}^{ag}=E_0^{ag}\left(1-D^{ag}\right)\\ E^{oitz}=E_0^{oitz}\left(1-D^{oitz}\right)\\ E^m=E_0^m\left(1-D^m\right)\\ E^{om}=E_0^{om}\left(1-D^{om}\right)\\ E^{\text{itz}}=E_0^{itz}\left(1-D^{itz}\right)\end{array}$ (6)

$D_t=\{\begin{array}{lll}0\hfill & \hfill & {\epsilon }_{\max }\le {\epsilon }_{t0}\hfill \\ 1-\frac{{\epsilon }_{t0}}{{\epsilon }_{\max }}+\frac{{\epsilon }_{\max }-{\epsilon }_{t0}}{{\eta }_t{\epsilon }_{t0}-{\epsilon }_{t0}}\frac{{\epsilon }_{t0}}{{\epsilon }_{\max }}\left(1-\mu \right)\hfill & \hfill & {\epsilon }_{t0}<{\epsilon }_{\max }\le {\eta }_t{\epsilon }_{t0}\hfill \\ 1-\frac{\mu }{{\xi }_t-{\eta }_t}\frac{{\epsilon }_{\max }-{\eta }_t{\epsilon }_{t0}}{{\epsilon }_{\max }}+\frac{\mu {\epsilon }_{t0}}{{\epsilon }_{\max }}\hfill & \hfill & {\eta }_t{\epsilon }_{t0}<{\epsilon }_{\max }\le {\xi }_t{\epsilon }_{t0}\hfill \\ 1\hfill & \hfill & {\epsilon }_{\max }>{\xi }_t{\epsilon }_{t0}\hfill \end{array}$ (7)

$D_c=\{\begin{array}{lll}1-\frac{\delta }{\omega }\hfill & \hfill & {\epsilon }_{\max }\le \lambda {\epsilon }_{c0}\hfill \\ 1-\frac{1-\delta }{1-\lambda }\frac{{\epsilon }_{\max }-\lambda {\epsilon }_{c0}}{{\epsilon }_{\max }}-\delta \frac{{\epsilon }_{c0}}{{\epsilon }_{\max }}\hfill & \hfill & \lambda {\epsilon }_{c0}<{\epsilon }_{\max }\le {\epsilon }_{c0}\hfill \\ 1-\frac{1-\omega }{1-{\eta }_c}\frac{{\epsilon }_{\max }-{\epsilon }_{c0}}{{\epsilon }_{\max }}-\frac{{\epsilon }_{c0}}{{\epsilon }_{\max }}\hfill & \hfill & {\epsilon }_{c0}<{\epsilon }_{\max }\le {\eta }_c{\epsilon }_{c0}\hfill \\ 1-\frac{\omega {\epsilon }_{c0}}{{\epsilon }_{\max }}\hfill & \hfill & {\eta }_c{\epsilon }_{c0}<{\epsilon }_{\max }\le {\xi }_c{\epsilon }_{c0}\hfill \\ 1\hfill & \hfill & {\epsilon }_{\max }>{\xi }_c{\epsilon }_{c0}\hfill \end{array}$ (8)

式中， ${\epsilon }_0$ 为峰值应变； $\eta$ 为残余应变系数； $\xi$ 为极限应变系数； $\lambda$ 为弹性应变系数； $\delta$ 为弹性抗压强度系数； $\omega$ 为残余抗压强度系数； $\mu$ 为残余抗拉强度系数；下标t和c分别代表单元抗拉、抗压两情况。

对于二维随机骨料模型，各相介质的体积分数可简化计算为面积分数。假设大网格划分尺寸为a，小网格划分尺寸为b(a > b)，大网格中有n个小网格单元属性判定为骨料，则骨料的面积分数c₁ = nb²/a²。同理分别用c₀，c₁，c₂，c₃，c₄表示新砂浆，骨料，老界面，老砂浆，新界面的面积分数。

基于等应力假设可得到

$\epsilon =\frac{\sigma }{E_m}=c_0\frac{\sigma }{E_0}+c_1\frac{\sigma }{E_1}+c_2\frac{\sigma }{E_2}+c_3\frac{\sigma }{E_3}+c_4\frac{\sigma }{E_4}=\sigma \left(\frac{c_0}{E_0}+\frac{c_1}{E_1}+\frac{c_2}{E_2}+\frac{c_3}{E_3}+\frac{c_4}{E_4}\right)$ (9)

串联等效单元的弹性模量为

$\frac{1}{E_{eq}}=\frac{c_0}{E_0^m\left(1-D^m\right)}+\frac{c_1}{E_0^{ag}\left(1-D^{ag}\right)}+\frac{c_2}{E_0^{oitz}\left(1-D^{oitz}\right)}+\frac{c_3}{E_0^{om}\left(1-D^{om}\right)}+\frac{c_4}{E_0^{itz}\left(1-D^{itz}\right)}$ (10)

对于再生混凝土，新老界面厚度较小，当采用均质化等效模型时网格较大，当网格大小大于老砂浆厚度，可能单一网格会包含多相介质，因此本文将再生混凝土其等效为三相：砂浆(m)、骨料(ag)、等效单元(em)，建立三相介质的细观等效化本构模型，其中 ${\epsilon }_0^{em}\le {\epsilon }_0^m\le {\epsilon }_0^{ag}\le {\epsilon }_r^{em}\le {\epsilon }_r^m\le {\epsilon }_r^{ag}\le {\epsilon }_u^{em}\le {\epsilon }_u^m\le {\epsilon }_u^{ag}$ 得到再生混凝土串联等效单元的弹性损伤模型的公式为

$\frac{1}{E^{eq}}=\{\begin{array}{l}\frac{c_0}{E_0^m}+\frac{c_1}{E_0^{ag}}+\frac{c_2}{E_0^{em}}\epsilon \le {\epsilon }_0^{em}\\ \frac{c_0}{E_0^m}+\frac{c_1}{E_0^{ag}}+\frac{c_2}{E_r^{em}}{\epsilon }_0^{em}<\epsilon \le {\epsilon }_0^m\\ \frac{c_0}{E_r^m}+\frac{c_1}{E_0^{ag}}+\frac{c_2}{E_r^{em}}{\epsilon }_0^m<\epsilon \le {\epsilon }_0^{ag}\\ \frac{c_0}{E_r^m}+\frac{c_1}{E_r^{ag}}+\frac{c_2}{E_r^{em}}{\epsilon }_0^{ag}<\epsilon \le {\epsilon }_r^{em}\\ \frac{c_0}{E_r^m}+\frac{c_1}{E_r^{ag}}+\frac{c_2}{E_u^{em}}{\epsilon }_r^{em}<\epsilon \le {\epsilon }_r^m\\ \frac{c_0}{E_u^m}+\frac{c_1}{E_r^{ag}}+\frac{c_2}{E_u^{em}}{\epsilon }_r^m<\epsilon \le {\epsilon }_r^{ag}\\ \frac{c_0}{E_u^m}+\frac{c_1}{E_u^{ag}}+\frac{c_2}{E_u^{em}}{\epsilon }_r^{ag}<\epsilon \le {\epsilon }_u^{em}\\ \frac{c_0}{E_r^m}+\frac{c_1}{E_r^{ag}}{\epsilon }_u^{em}<\epsilon \le {\epsilon }_u^m\\ \frac{c_1}{E_r^{ag}}{\epsilon }_u^m<\epsilon \le {\epsilon }_u^{ag}\\ 0{\epsilon }_u^{ag}<\epsilon \end{array}$ (11)

式中 ${\epsilon }_0$ 、 ${\epsilon }_r$ 、 ${\epsilon }_u$ 分别为峰值应变、残余应变和极限应变。 ${\epsilon }_r=\eta {\epsilon }_0$， ${\epsilon }_u=\xi {\epsilon }_0$。

本文考虑了再生混凝土各相中材料的随机分布，各参数取值见表1。

Table 1. Material parameter value

表1. 材料参数取值

3. 再生混凝土立方体单轴静态压缩数值模拟

单轴压缩加载模型

本文运用基于余能原理的基面力元法 [4]，编制考虑均质化模型的计算分析程序，建立立方体再生混凝土试件加载模型，试件选取100 mm × 100 mm × 100 mm，进行单轴压缩数值模拟试验，首先将立方体模型简化为二维模型取截面尺寸为100 mm × 100 mm，加载模型如图3。加载时采用竖向加载，加载采用逐级静力位移加载，每级加载位移为0.01 mm。

Figure 3. Load model

图3. 加载模型

随机骨料模型选取不同的随机数进行投放骨料可以得到骨料分布位置不同的试件，本文选取三个由三组骨料分布不同、骨料颗粒数相同的二维数值模型，试件如图4。

Figure 4. Two-dimensional diagram of homogenized model

图4. 均质化模型试件二维图

运用余能原理基面力元法细观损伤计算程序计算分析，对生成的三个随机试件的串联等效的均质化模型进行单轴压缩试验数值模拟。三个试件计算结果和试验数据 [15] 列于表2中。同时以应变为横坐标，应力为纵坐标绘制应力–应变全曲线，串联等效的计算结果如图5。

由图5可以看出，再生混凝土的均质化模型压缩加载过程中再生混凝土试件在初始阶段处于弹性阶段，应力随着压缩的应变增大而增大，当应力趋近于峰值时开始上升缓慢。峰值过后，进入损伤阶段，应力随应变的增加而降低，直至稳定于一个较低的值，不再变化。这与再生混凝土各相本构模型的抗压部分的规律一致。

Table 2. Numerical simulation results data

表2. 数值模拟计算结果数据

Figure 5. Uniaxial compression stress-strain curve of series model

图5. 串联模型单轴压缩应力–应变曲线

为明显简洁的得到再生混凝土数值模拟计算的损伤破坏图，采用FORTRAN中的Quick Win模块将各单元的不同阶段以不同的颜色显示出来。等效单元将设为同一种颜色黄绿色，砂浆为橘色，骨料为蓝色，破坏单元用黑色块来表示，观察数值模拟试件的破坏形态，如图6。

由损伤图分析试件在单轴压缩加载过程中损伤状态的变化，可以得出再生混凝土的破坏裂纹规律：加载到一定应变时，试件局部单元进入损伤阶段，然后沿一定角度逐渐展开，直至整个试件破坏。破坏形态与普通混凝土相近，裂纹相连最后形成四角锥形状；再生混凝土压缩试验 [15] 中破坏时斜裂缝与加载方向成58˚~64˚角，这与数值模拟分析得到的损伤图一致。由损伤图可以看出，损伤开始出现的位置一般为再生骨料周围，这是因为骨料周围的等效单元包含强度较低的老界面和新界面，最易先达到损伤阶段。均质化等效模型试件破坏裂缝的方向与发展与随机骨料模型一致，并与实际试验规律一致。

Figure 6. Uniaxial compression numerical simulation damage diagram

图6. 单轴压缩数值模拟计算损伤图

4. 结论

1) 利用均质化分析方法，推导出均质化模型的串联等效本构关系，建立了再生混凝土材料的多折线损伤本构模型。

2) 基于余能原理基面力元法，开发了再生混凝土均质化分析的非线性基面力元分析软件和均质化前处理软件。

3) 采用串联等效的均质化模型对再生混凝土的单轴压缩试验进行了数值计算分析，得到了应力–应变软化曲线和损伤破坏过程。

4) 通过将本方法的结果与试验数据的对比，验证了模型建立的可行性与合理性。

5) 研究表明，采用均质化模型的计算效率大大提高，远远高于基于随机骨料模型的细观损伤分析方法，且能保证一定的计算精度。