基于隐式梯度理论的混凝土细观随机断裂模型

摘要：根据多尺度力学的基本思想，将隐式梯度非局部化理论与串并联弹簧模型理论相结合，建立了混凝土细观随机断裂模型。通过引入微弹簧断裂应变为服从对数正态分布的随机变量，验证了混凝土材料微观断裂的不可控制性；通过引入非局部化比例，解决了传统串并联弹簧模型的层数敏感性问题。对建立的模型从特征参数、局部应变与损伤的演化发展和层数敏感性几个方面进行了分析。最后，将模型结果与试验结果进行了对比，验证了该模型的正确性。研究结果表明：该模型计算效率高，且更加有利于在计算软件中实现。

关键词：混凝土；隐式梯度；非局部化；细观随机断裂模型

中图分类号：TU528.01文献标志码：A

0引言

从材料组成成分看，混凝土材料是多组分、多相和非均质的复合材料，其内部组织结构包含水泥石、不同形状和大小的集料以及各种毛细管空隙结构等，导致混凝土的宏观力学指标具有显著随机性[1]。同时，在单轴受力情况下，混凝土表现出应变软化、强度（刚度）退化、塑性变形、残余应力等特征，应力应变全过程曲线具有明显的非线性特点[2]。如何反映混凝土受力力学行为的非线性与随机性，构成了混凝土力学研究的核心科学问题。

1990年，Krajcinovic等[3]在此基础上假设弹簧破坏强度为服从相同分布的随机变量，采用概率定义损伤建立了混凝土单轴受力的细观弹簧模型。由于Krajcinovic等定义的损伤是基于细观弹簧破坏概率，因此该模型给出的是确定性的损伤变量。1996年，Kandarpa等[4]在Krajcinovic模型基础上做出改进，通过随机场的相关结构考虑了细观弹簧之间的相互影响，建立了基于细观弹簧模型的混凝土随机损伤模型。但是该模型存在2个方面的问题：其一，以受拉破坏为理论依据却采用受压试验结果建模；其二，在模型中以变形代替了应变，并不具有真正的细观意义。1999年以来，李杰等[1，5]在细观层次上将混凝土离散为具有一定特征高度和截面积的小柱体，并用微弹簧模型表示，且假定微弹簧的断裂应变服从某一随机场，从而建立了细观随机断裂模型。较之传统模型，这一类模型能从细观上统一解释混凝土应力应变曲线的随机性和非线性特征，实现了随机性与非线性的综合反映，为多维随机损伤本构关系[6]的建立奠定了基础。然而，进一步研究表明，在采用细观随机断裂模型通过串并联弹簧方式模拟混凝土单轴受力应力应变关系时，模型的计算结果对微弹簧层数的选取有很大的依赖性。这一问题的本质，在于经典本构关系在涉及软化问题时所表现出的局部性限制。为解决这一问题，2012年，刘汉昆等[7]将积分型非局部化理论与上述串并联弹簧相结合，引入积分型非局部化模型，对控制损伤的变量进行非局部化处理，通过将控制损伤的变量在一定范围内进行加权平均，从而考虑损伤的非均匀性[8]。在有限元实现过程中，这类模型要求光滑长度必须大于网格尺寸，非局部变量才可以求解，这实际上是高估了微结构间的相互作用，而且直接积分形式更多的时候只能直接求得一些特定点的非局部应变，要求解任意的非局部变量需要专门的数值技巧[9]，因此积分型非局部化模型应用于实际结构的数值模拟存在不足。

有鉴于此，本文中引入梯度型非局部化理论，基于多尺度力学的基本思想，在细观层次引入断裂的随机描述，建立了基于梯度型非局部化理论的混凝土随机断裂本构模型。较之Kandarpa等建议的模型，本文模型能真实客观地反映混凝土随机损伤应力重分布损伤演化的物理过程，同时，该模型中的每一个微弹簧取为细观单元，其应力应变关系为非线性关系，不再是线性关系[8]。较之积分型非局部化模型，该模型的非局部化过程只依赖于所考虑材料点相关变量的梯度，不再需要计算整个结构区域的所有单元信息，计算效率得到了很大的提高，且更加有利于在计算软件中实现。

1细观弹簧模型

为了模拟混凝土单轴受力试件的力学行为，采用如图1所示的串并联弹簧模型，其中，N×M个基本单元中的每一个都用一微弹簧束建模，M个微弹簧束两端与刚性板相连接形成层单元，N个弹簧层单元体串联形成整个体系。图1中F为轴力。

各个微弹簧的断裂应变为随机变量，故可以反映混凝土材料微观断裂的随机性。

当外力小于某一值时，微弹簧系统不发生断裂，相应的细观应力应变关系具有线性变化特征。当外力增加，使得具有最小断裂应变的微弹簧的应变达到其断裂界限值时，该微弹簧断裂，所释放的应力由其他微弹簧承担，形成应力重分布，所以细观随机断裂模型在物理上科学反映了混凝土受力力学行为的随机性和非线性。

根据上述模型，第i层微弹簧束（细观层次）损伤变量Di可定义为

式中：Ai为第i个细观单元的横截面积；ADi为第i个细观单元中微弹簧断裂而导致材料退出工作的面积。

假定材料离散后模型中的细观单元横截面积均相等，损伤变量Di可定义为

式中：H（x）为Heaviside函数；Aji为第i个细观单元中第j个微弹簧面积；A为宏观单元的横截面积，A=Ai；εjei为第i个细观单元中第j个微弹簧的弹性应变；Δji为第i个细观单元中第j个微弹簧的极限破坏应变，为服从对数正态分布的随机变量。

当微弹簧个数M+∞时，并联微弹簧可以看作一维连续体，式（2）的极限Di（εe）为

式（1）～（4）的建立运用了2个尺度的物理关系，即微观层次的随机断裂、细观层次的应力应变关系。而对于宏观尺度，根据力的平衡和位移协调条件可得

式中：σMA，εMA分别为串并联弹簧体系的宏观平均应力和宏观平均应变；E0为材料弹性模量；σi为第i层弹簧的总应力；εi为第i层弹簧的总应变；εpi为第i层弹簧的塑性应变。

对于单轴受拉模型，取εpi=0；对于单轴受压模型，εpi可以按经验公式进行取值，即塑性应变εpi=ξp[Di（εei）]pεei[10]，p为塑性参数，ξp为塑性应变参数。

值得指出的是，按照细观随机断裂模型分析宏观单轴受力试件，无论是单轴受拉还是单轴受压，均存在层数敏感性问题[1114]。图3为根据上述模型计算的1个单轴受拉试件的试验结果，其中各参数取值为：材料弹性模量E0=42 GPa；断裂应变随机场参数的均值λ=4.92；标准差ζ=0.80；塑性应变参数ξp=0；塑性参数p=0（对应于单轴受拉工况）。

图3层数敏感性

Fig.3Layer Sensitivity从图3可以看出，计算结果严重依赖于模型所选取的弹簧层数，弹簧层数N越大，计算结果的下降段越陡峭。

图4为串并联弹簧模型的简化模型，其中，ft为抗拉强度，ε0为抗拉强度对应的峰值应变，Et为线性软化段斜率。现以图4所示的简化模型做进一步分析，可以更直观地解释上述层数敏感性问题，其中每个“圆球”代表一个细观单元，并假定图2（a）中宏观单元的总长度为L0（共有n个细观单元）。为方便后续分析，将图2（b）中细观单元的本构模型简化为线性软化函数。

由于初始缺陷的影响，某个单元开裂后，应变局部化导致变形集中于断裂单元，而其他n-1个未断裂单元则为线弹性，此时

图5为数值分析结果。从图5可以看出，采用局部应力应变关系计算得到的应力应变结果与所采用的单元数量密切相关。特别当n+∞时，整体的应力应变曲线退化为线弹性，这意味着整个破坏过程中的能量耗散为0，显然这与其物理过程的事实不符合。

图5数值分析结果

Fig.5Simulating Analytical Results上述问题的根源在于细观本构关系在本质上属于局部本构关系。而对于混凝土等软化材料而言，大量存在的微裂纹之间会产生相互作用，使得材料变形呈现出高度的应变局部化现象，从而使局部本构模型的独立性假定难以满足，相应的有限元分析结果会严重依赖于网格尺寸，即所谓的网格敏感性问题[1516]。3混凝土非局部化损伤模型

为了能合理地反映混凝土材料的应变软化现象并解决应变局部化引起的数值问题，建立非局部化本构模型是较为理想的选择。非局部化本构模型分为积分型非局部化本构模型和梯度型非局部化本构模型[9]。

3.1积分型非局部化本构模型

积分型非局部化本构模型[1718]的基本思想是将与材料本构关系有关的物理量在一定空间内进行加权平均，从而考虑材料细观结构的非均匀性，则有

（x）=∫1Vα（x，ξ）f（ξ）dξ（13）

式中：（x）为非局部变量；f（ξ）为对应的局部变量；V为局部平均单元体积；α（x，ξ）为平均化权函数，满足归一化条件∫1Vα（x，ξ）dξ=1。

3.2梯度型非局部化本构模型

将上述积分型非局部化本构模型进行Taylor展开，可以将变量的高阶梯度引入材料本构关系中，给出一类梯度型非局部模型[1920]。

将局部变量在x位置进行Taylor展开，即

f（ξ）=f（x）+f（x）（ξ-x）+1122f（x）・

（ξ-x）2+1132f（x）（ξ-x）2+…（14）

将式（14）代入式（13），可以得到

（x）=f（x）+c2f（x）+d4f（x）（15）

式中：c，d均为参数，与相关函数α（x，ξ）有关；f（x）为局部变量。

若仅保留二阶梯度项，可得到如下常见的显式非局部梯度模型

（x）=f（x）+c2f（x）（16）

由式（16）可见，在显式非局部梯度模型中，非局部变量（x）由其对应的局部变量f（x）及其二阶梯度2f（x）显式地确定。

上述显式非局部梯度模型要求给出局部变量的二阶梯度，其有限元插值格式的构造需要高阶连续插值函数，这给模型的有限元实现带来了困难。为此，一般将其转化为如下隐式梯度模型

（x）-c2（x）=f（x）（17）

在求解上述非局部化隐式梯度模型时，必须补充新的边界条件。一般可采用如下自然边界条件

（x）T=0（18）

式中：T为计算域边界的外法线矢量。

研究结果表明[17]，对于损伤模型或软化塑性模型等非线性模型，非局部变量只能选择在整个加载过程中不会发生减小的内变量，如损伤变量、损伤能释放率、累积塑性应变等，故本文中选取局部弹性应变为非局部变量（对于反复加载情况，这里的弹性应变应取为历史加载最大值）。

3.3数值实现

本文中在细观单元尺度构造有限元方法求解隐式梯度方程。由有限元基本知识，可以得到式（16）的等效积分弱形式，即

∫1LW（x）dx-∫1LWc2（x）dx=

∫1LWf（x）dx（19）

式中：L为归一化特征长度。

根据分部积分，则有

∫1LcW2（x）dx=cW（x）|L2L1-

c∫1LW（x）dx（20）

根据边界条件cW（x）|L2L1=0，将式（20）代入式（19），则有

∫1LW（x）dx+c∫1LW（x）dx=

∫1LWf（x）dx（21）

按照有限元插值理论，对节点进行插值，即

（x）=J

f（x）=Jg

W=Jw（22）

式中：J为形函数；w为参数。

将式（22）代入式（21），可得

（∫1LJTJdx+c∫1LBTBdx）=（∫1LJTJdx）g（23）

式中：B=J1x。

取形函数J为三节点二阶的拉格朗日插值函数，则有

Ji（x）=（x-xi-1）（x-xi+1）1（xi-xi-1）（xi-xi+1）（24）

J=（J1，J2，J3）

g=（g1，g2，g3）（25）

将式（24）代入式（23）可得

A=Bg（26）

式中：A，B均为3×3阶的系数矩阵。

系数矩阵A，B的具体表达式为

A（i，j）=∫1L（Ji（x）Jj（x）+cBi（x）Bj（x））dx（27）

B（i，j）=∫1L（J21（x）+cB21（x））dxi=j=1

∫1LJi（x）Jj（x）dx其他（28）

至此，非局部化过程的数值实现已经转化为线性方程式（26）的求解。

事实上，拉格朗日插值函数并不是阶次越高精度就越高，高阶插值可能引起数值震荡，即龙格现象，因此，本文中采用分段二阶拉格朗日插值。4非局部化模型分析

若非注明，以下分析采用的参数为：材料弹性模量E0=42 GPa，断裂应变的均值λ=4.92；标准差ζ=0.80；塑性应变参数ξp=0；p=0（对应单轴受拉工况）。

4.1参数c的讨论

文献[21]中指出，参数c与材料特征长度l之间存在以下关系

c=l2/2（29）

式中：l取3倍骨料粒径[18]。

图6为参数c=0.001，0.01，0.05，0.1，0.2时对应的应力应变（σε）曲线。

图6不同c值对应的σε曲线

Fig.6σε Curves Corresponding to Different c从图6可以看出，当c=0.001时，应力跌落现象十分明显，此时存在严重的应变局部化现象。随着c值的增大，对应的非局部化程度也会增大，应力跌落现象能得到明显改善。

4.2应变局部与损伤演化发展

为了分析不同c值对应的细观应变与损伤演化，分别取c=0.000 8，0.001，0.001 2，考察宏观应变ε=1.5×10-4时细观应变和细观损伤沿宏观方向的分布，如图7所示，其中，D为损伤指数。

图7ε=1.5×10-4时应变分布与损伤演化

Fig.7Strain Distributions and Damage

Evolution when ε=1.5×10-4从图7可以看出，非局部化损伤模型中的损伤和应变集中在一定宽度的窄带内，随着c值的增大，细观应变和损伤集中的范围也逐渐增大，对应的非局部化程度越高，应变局部化程度削弱。

图8为c=0.01时宏观应变分别为ε1=1.5×10-4，ε2=2×10-4，ε3=2.3×10-4，ε4=2.7×10-4，ε5=3×10-4，ε6=3.3×10-4的细观应变发展和损伤演化曲线。

图8c=0.01时细观应变与损伤演化

Fig.8Meso Strain Distributions and Damage

Evolution when c=0.01从图8可以看出，当宏观应变较小时，材料尚未进入软化段，此时应变和损伤都是均匀发展的。随着宏观应变的增大，材料进入软化段，此时由于软化段应变局部化的影响，应变和损伤沿宏观长度方向的分布将不再均匀，而且宏观应变越大，不均匀程度越明显。

4.3层数敏感性

图9中给出了c取不同的值，弹簧层数N取4个不同值对应图9不同c值对应的层数敏感性计算结果

Fig.9Calculation Results of Layer Sensitivety

Corresponding to Different c Values的计算结果。结果表明，代表非局部化程度的参数c一旦确定，模型得到的应力应变曲线将不依赖于模型层数，不同层数对应的计算结果相差甚小，这说明本模型克服了层数敏感性问题。

5算例分析

为了验证本文所建立的模型，分别对单轴受拉和单轴受压试验进行模拟，其中，单轴受拉试验数据来源于文献[10]，单轴受压试验数据来源于笔者重新做的一批C30混凝土棱柱体试验，试件尺寸为150 mm×150 mm×300 mm。根据试验数据，用参数识别的方式得到微观断裂应变的随机场参数（单轴受压工况还需要识别塑性应变参数），各参数取值如表1所示。图10，11分别为模拟结果与单轴受拉、单轴受压试验结果对比。

表1各参数取值

Tab.1Values of Each Parameter参数1E0/GPa1λ1ζ1ξp1p1c单轴受拉138.2014.9610.5210.0010.0010.08单轴受压135.8017.2310.6610.5310.4910.20图10单轴受拉试验结果与模拟结果的对比

Fig.10Comparison of Test Result and Simulation

Result of Uniaxial Tension图11单轴受压试验结果与模拟结果的对比

Fig.11Comparison of Test Result and Simulation

Result of Uniaxial Compression从模拟结果与试验结果的对比可以看出，本文中建议的模型能较好地模拟混凝土单轴受拉、单轴受压试验。相对单轴受压试验，单轴受拉试验参数c的取值较小，对应的非局部化程度较小，应力跌落现象较明显。6结语

以多尺度分析的基本思想为基础，建立了具有微观、细观、宏观3个尺度的混凝土随机损伤模型，为反映混凝土应力应变曲线的随机性和非线性提供了物理背景。在此基础上，本文中引入隐式梯度非局部化理论，考虑了混凝土材料在变形后期存在的变形局部化问题，解决了局部化弹簧模型中存在的层数敏感性问题。最后通过模拟结果与试验结果的对比可以看出，本文中建立的模型是合理的。参考文献：

References：[1]李杰，卢朝辉，张其云.混凝土随机损伤本构关系――单轴受压分析[J].同济大学学报：自然科学版，2003，31（5）：505509.

LI Jie，LU Zhaohui，ZHANG Qiyun.Study on Stochastic Damage Constitutive Law for Concrete Material Subjected to Uniaxial Compressive Stress[J].Journal of Tongji University：Natural Science，2003，31（5）：505509.

[2]李杰.混凝土随机损伤力学的初步研究[J].同济大学学报：自然科学版，2004，32（10）：12701277.

LI Jie.Research on the Stochastic Damage Mechanics for Concrete Materials and Structures[J].Journal of Tongji University：Natural Science，2004，32（10）：12701277.

[3]KRAJCINOVIC D，STOJIMIROVIC A.Deformation Processes in Semibrittle Polycrystalline Ceramics[J].International Journal of Fracture，1990，42（1）：7386.

[4]KANDARPA S，KIRKNER D J.Stochastic Damage Model for Brittle Materiel Subjected to Monotonic Loading[J].Journal of Engineering Mechanics，1996，122（8）：788795.

[5]李杰，张其云.混凝土随机损伤本构关系[J].同济大学学报：自然科学版，2001，29（10）：11351141.

LI Jie，ZHANG Qiyun.Study of Stochastic Damage Constitutive Relationship for Concrete Material[J].Journal of Tongji University：Natural Science，2001，29（10）：11351141.

[6]LI J，REN X.Stochastic Damage Model for Concrete Based on Energy Equivalent Strain[J].International Journal of Solids and Structures，2009，46（11）：24072419.

[7]刘汉昆，李杰.单轴受拉状态下混凝土非局部化细观损伤模型[J].同济大学学报：自然科学版，2014，42（2）：203209.

LIU Hankun，LI Jie.Nonlocal Microdamage Model for Concrete Under Uniaxial Tension[J].Journal of Tongji University：Natural Science，2014，42（2）：203209.

[8]李杰，吴建营，陈建兵.混凝土随机损伤力学[M].北京：科学出版社，2013.

LI Jie，WU Jianying，CHEN Jianbing.Stochastic Damage Mechanics of Concrete Structures[M].Beijing：Science Press，2013.

[9]汪忠明.基于非局部增强梯度模型的混凝土断裂损伤研究[D].合肥：合肥工业大学，2010.

WANG Zhongming.Research on Fracture and Damage Concrete Based on the Enhanced Damage Model[D].Hefei：Hefei University of Technology，2010.

[10]曾莎洁.混凝土随机损伤本构模型与试验研究[D].上海：同济大学，2012.

ZENG Shajie.Experiment Research on Model for Concrete[D].Shanghai：Tongji University，2012.

[11]李朝红，王海龙，徐光兴.混凝土梁受弯破坏的三维细观数值试验[J].中国公路学报，2010，23（2）：5055.

LI Chaohong，WANG Hailong，XU Guangxing.Mesoscopic 3D Numerical Test of Bending Fracture of Concrete Beam[J].China Journal of Highway and Transport，2010，23（2）：5055.

[12]李建波，林皋，陈健云，等.混凝土损伤演化的随机力学参数细观数值影响分析[J].建筑科学与工程学报，2007，24（3）：712.

LI Jianbo，LIN Gao，CHEN Jianyun，et al.Mesoscopic Numerical Influence Analysis of Random Mechanical Parameter of Concrete Damage Evolution[J].Journal of Architecture and Civil Engineering，2007，24（3）：712.

[13]林皋，李建波，赵娟，等.单轴拉压状态下混凝土破坏的细观数值演化分析[J].建筑科学与工程学报，2007，24（1）：16.

LIN Gao，LI Jianbo，ZHAO Juan，et al.Mesoscopic Numerical Evolution Analysis of Concrete Damage Under Uniaxial Tension and Compression[J].Journal of Architecture and Civil Engineering，2007，24（1）：16.

[14]丰茂东，李建波，林皋，等.随机力学参量对混凝土细观损伤演化的影响[J].建筑科学与工程学报，2010，27（4）：16.

FENG Maodong，LI Jianbo，LIN Gao，et al.Influence of Stochastic Mechanical Parameters on Mesodamage Evolution of Concrete[J].Journal of Architecture and Civil Engineering，2010，27（4）：16.

[15]RUDNICKI J W，RICE J R.Conditions for the Localization of Deformation in Pressuresensitive Dilatant Materials[J].Journal of the Mechanics and Physics of Solids，1975，23（6）：371394.

[16]FEENSTRA P H，DE BORST R，ROTS J G.Numerical Study on Crack Dilatancy.Part 2：Applications[J].Journal of Engineering Mechanics，1991，117（4）：754769.

[17]BAZANT Z P，PIJAUDIERCABOT G.Nonlocal Continuum Damage，Localization Instability and Convergence[J].Journal of Applied Mechanics，1988，55（2）：287293.

[18]BAZANT Z P，JIRASEK M.Nonlocal Integral Formulations of Plasticity and Damage：Survey of Progress[J].Journal of Engineering Mechanics，2002，128（11）：11191149.

[19]PEERLINGS R H J，DE BORST R，BREKELMANS W A M，et al.Gradient Enhanced Damage for Quasibrittle Materials[J].International Journal for Numerical Methods in Engineering，1996，39（19）：33913403.

[20]PEERLINGS R H J，DE BORST R，BREKELMANS W A M，et al.Gradientenhanced Damage Modelling of Concrete Fracture[J].Mechanics of Cohesivefrictional Materials，1998，3（4）：323342.

[21]PEERLINGS R H J，GEERS M G D，DE BORST R，et al.A Critical Comparison of Nonlocal and Gradientenhanced Softening Continua[J].International Journal of Solids and Structures，2001，38（44）：77237746.《工业建筑》2015年征订通知

《工业建筑》由中国钢铁工业协会主管，中冶建筑研究总院有限公司主办，1964年创刊；重点报道建筑设计、建筑结构、地基和基础、建筑材料、施工技术等方面的应用研究成果及开发技术，突出工业建筑、钢结构、鉴定加固技术等特点；优先报道国家或者省部科学基金资助项目和重大工程建设项目；内容丰富翔实，理论结合实际；适合土木工程界广大科研、设计、教学、施工等专业技术人员及高校师生阅读。

《工业建筑》为大16开本，128页，每月20日出版，国内外公开发行，每期定价15元，全年共180元。全国各地邮局均可订阅，邮发代号：2825，也可直接汇款到本编辑部订阅。

《工业建筑》兼营广告，每年8月开始征订下一年广告，收费合理，具有一定规模，是展示企业形象的理想舞台。广告许可证：京海工商广字8231。

地址：北京市海淀区西土城路33号邮编：100088

Email：.cn网址：http：//

电话：（010）82227677（广告），（010）82227237/7239/8043/8927（编辑），

（010）82227236（发行），（010）82227008/7239（传真）