裂缝宽度验算是混凝土结构正常使用极限状态验算的重要内容之一。对于一些受力条件单一、边界条件明确的简单钢筋混凝土构件(譬如：钢筋混凝土轴心受拉构件、钢筋混凝土受弯构件等)，过去的许多研究以及规范都建议了半经验半理论的裂缝宽度简化计算公式[1]。然而，在实际结构体系中，混凝土结构构件的内力状态和边界条件都比较复杂，构件类型也十分丰富(尤其是近年来，许多新型钢-混凝土组合构件以及新材料-混凝土组合构件不断涌现)，亟需发展更为通用的混凝土裂缝宽度计算方法以适应各类混凝土结构构件正常使用极限状态的验算需求。
有限元方法是目前结构工程设计和研究中应用最普遍的数值方法，采用有限元方法计算结构构件中混凝土的裂缝宽度不仅通用性强，也易于被结构工程师所接受。经典的弥散裂缝模型(Smeared Crack Model)以及Bazant和Oh[2]于1983年提出的裂缝带理论(Crack Band Theory)为采用有限元方法进行混凝土裂缝分析提供了重要的理论基础，至今，已有大量学者基于这些理论背景采用基于弥散裂缝模型的通用有限元程序(例如：ABAQUS、MSC.MARC、ANASYS等)开展了一系列丰富的混凝土裂缝分析的实践，然而，这些实践几乎均以“开裂应变”作为最终结果的展示[3—10]，虽然能够定性地描述混凝土结构构件的开裂程度以及裂缝分布情况，但却无法定量地给出裂缝宽度值，不仅使计算结果的准确性难以得到定量的检验，也无法满足结构设计规范中“裂缝宽度验算”的需求。此外，在分析中如何合理地进行参数设置至今也未有成熟的研究结论。分层壳单元是一种近年来发展迅速的适用于钢筋混凝土板、剪力墙等构件精细化模拟的高效数值模型，在工程实践中应用广泛。然而，如何应用分层壳模型实现钢筋混凝土构件裂缝宽度的定量预测，仍鲜有研究。
在上述背景下，本文的主要目标是：1) 以弥散裂缝模型和裂缝带理论为基础，给出采用通用有限元程序的分层壳单元计算配筋混凝土构件裂缝宽度的计算流程；2) 以某承受负弯矩的钢-混凝土简支组合梁为算例，识别影响混凝土裂缝宽度计算结果的关键参数，从而对计算过程中的合理参数设置提出建议。需要特别说明的是，为便于实际工程的应用，本文的讨论均不涉及自编程序或基于通用有限元程序的二次开发，而是讨论如何采用已有成熟的通用有限元程序中的弥散裂缝模型和分层壳模型合理地完成配筋混凝土构件的裂缝宽度计算。所有算例均采用 MSC.MARC[11]完成，但相关方法和结论对于当前工程和科研中广泛采用的其他集成弥散裂缝模型的通用有限元程序(例如：ABAQUS和ANASYS等)均适用。
1 计算理论与流程
1.1 基本假定
本文采用以下三个基本假定：
第一，本文的研究涉及“裂缝宽度”这一非力学物理量和网格尺寸之间的讨论，所有这些讨论的前提是：网格密度首先要满足位移、反力、应力、应变等力学物理量的求解精度，不可过于稀疏。
第二，本文仅研究分层壳单元，不直接模拟钢筋与混凝土之间的粘结滑移行为，钢筋与混凝土完全变形协调，这也是钢筋混凝土结构分析中常用的假定。其他建模方式不在本文讨论范围内。
第三，本文讨论的通用有限元程序需集成弥散裂缝模型框架，对于没有集成该模型框架的通用有限元程序，无法采用本文建议的方法。
1.2 基于弥散裂缝模型的经典裂缝带理论概述及其局限性
在Bazant和Oh[2]所提出的经典裂缝带理论中，最基本的假定就是：采用断裂能Gf来描述一根裂缝开裂全过程的力学行为。断裂能Gf的定义是：一根裂缝开裂全过程单位面积所耗散的能量，它代表混凝土材料单轴拉应力σ与裂缝宽度w关系曲线的下包面积。之所以引入断裂能Gf的概念，是因为大量的混凝土直拉试验表明，同一种混凝土材料无论直拉试件的尺寸如何变化，其测得的σ-w曲线下包面积，也就是断裂能 Gf，基本是恒定的，这就充分证明了断裂能Gf是一种混凝土材料性能参数，与其他诸如立方体抗压强度 fcu、棱柱体抗压强度 fc、圆柱体抗压强度fc′等常用的混凝土材料性能参数具有同等的位置。笔者团队[12]曾采用文献[13—17]中的 17个不同参数的试件对比了Bazant和Oh[2]、Wittmann等[13]、欧洲模式规范CEB-FIP MC90[18]建议的断裂能 Gf公式的计算效果，结果表明欧洲模式规范CEB-FIP MC90[18]公式的计算结果与试验结果吻合最好，该公式建议断裂能Gf主要与混凝土圆柱体抗压强度 fc′以及最大骨料粒径 Dmax有关，具体形式为：
式中： fc′的单位是 N/mm2，若采用中国规范的混凝土棱柱体抗压强度 fc时，应进行强度转换 fc =0.95fc′；Gf的单位是N/mm；α为系数，α与最大骨料粒径 Dmax有关，Dmax=8 mm 时，α=0.025；Dmax=16 mm时，α=0.03；Dmax=32 mm 时，α=0.058。
图1所示为裂缝带理论的一些重要概念。当确定了混凝土轴心抗拉强度ft和断裂能Gf后，只需假定σ -w的曲线形式(通用有限元程序中的弥散裂缝模型常采用最简单的线性软化形式，如图 1(a)所示)，就可以唯一地确定混凝土材料的σ -w关系。然而，有限元计算是以材料的应力-应变关系(σ -ε)为基础，而σ -w关系虽然能够贴切地刻画混凝土材料的开裂特性，却无法直接用于有限元计算，因此将σ -w关系转换为σ -ε关系则成为非常关键的一步，而这一转换的核心就是将一根集中的裂缝宽度w在一个特定的范围内均匀弥散成应变ε，这就是“弥散裂缝”的概念，而这个特定的弥散范围就是所谓的“裂缝带”。因此，经典裂缝带理论中的裂缝带宽hc其实就是裂缝宽度向应变弥散的标距：
式中，εcr为开裂应变，代表一根裂缝宽度在裂缝带范围内的弥散化。
有了式(2)这一转换关系，σ -w关系就可以转换为用于有限元计算的σ -ε关系，当采用如图 1(b)所示最简单的受拉线性软化形式时，Bazant和Oh[2]最终给出了受拉开裂后软化模量 Ets的计算公式如式(3)所示，在这个公式中，除了裂缝带宽hc，其他所有参数都是材料属性参数。
式中：ft为混凝土轴心抗拉强度；Ec为混凝土弹性模量。
图1 裂缝带理论的一些重要概念
Fig.1 Some important concepts in the crack band theory
当完成有限元分析后，得到的开裂应变εcr还需按照式(2)乘以裂缝带宽 hc才能得到最终需要的裂缝宽度w值。
综上所述，式(2)和式(3)是裂缝带理论的两个核心方程，在这两个方程中，裂缝带宽hc无疑是最关键的参数。因此，采用裂缝带理论进行混凝土开裂有限元分析的核心就是要选取合适的裂缝带宽hc，经典的裂缝带理论认为裂缝带宽 hc应选为单元特征尺寸lele(垂直于裂缝方向单元尺寸的投影)，也就是当采用不同的单元网格尺寸分析同一个问题时，如果采用相同的参数设置，那么计算结果就与网格相关，如果要使计算结果一致，那么参数取值就要与网格相关。
网格相关是经典裂缝带理论最重要的结论之一，其背后的力学本质是混凝土开裂作为一种软化行为具有显著的局部化(Localization)特点。为了更清晰地阐明这一机理，本文采用通用有限元程序MSC.MARC[11]中的弥散裂缝模型完成了一个算例。图2所示为算例的基本信息及计算结果。
该算例为跨中承受竖向集中荷载的一根简支无筋素混凝土梁，该梁的基本尺寸如图2(a)所示，采用分层壳单元模拟该梁的全过程开裂行为。分别建立三个数值模型：模型A采用100 mm的单元网格进行计算，计算过程中的裂缝带宽hc也取为100 mm，与单元尺寸相同；模型B在模型A的基础上加密1倍网格，采用50 mm的单元网格进行计算，但计算过程中的裂缝带宽hc仍保持不变，取为100 mm；模型C同样在模型A的基础上加密一倍网格，采用50 mm的单元网格进行计算，但计算过程中的裂缝带宽hc也随之变为与单元网格50 mm一致。表1详细列出了这三个模型的具体参数取值，混凝土材料受压性能设为弹性，只考虑其受拉开裂的行为。采用位移控制对模型进行加载，分为均匀的400个加载子步逐步对跨中施加竖向位移至2 mm。采用残余力的收敛准则，收敛误差限设为0.02。
图2 跨中单点加载素混凝土梁算例Fig.2 Example of plain concrete beam subjected to mid-span point load
表1 素混凝土梁数值算例参数
Table 1 Parameters of numerical example of plain concrete beam
模型号单元网格尺寸lele/mm无缝带宽hc/mm软化模量Ets/(N/mm2)断裂能Gf/(N/mm)混凝土轴心抗拉强度ft/(N/mm2)混凝土弹性模量Ec/(N/mm2)材料泊松比νc剪力传递系数η模型A 100 100 3000模型B 50 100 3000 0.165 3 30000 0.2 0.1模型C 50 50 1428.6
图2还给出了三个模型的计算结果对比，包括：跨中梁底最大裂缝宽度w随跨中竖向位移δ 的变化曲线(如图2(b)所示)以及跨中荷载P随跨中梁底最大裂缝宽度w的变化曲线(如图2(c)所示)，在计算单元裂缝宽度w的过程中，提取的开裂应变为单元4个积分点开裂应变值的平均值。这些计算结果可从以下两个方面进行讨论：
1) 对比模型 A和模型 B，两者采用的单元网格尺寸不同，而计算参数hc取值相同，所得的计算结果差距很大，由此可证明：当采用不同的单元网格尺寸分析同一个问题，如果采用相同的参数设置，那么计算结果就与网格相关。
2) 对比模型 A和模型 C，两者采用的单元网格尺寸不同，而计算参数hc也随着网格尺寸的变化按照经典裂缝带理论进行调整，所得的计算结果几乎相同，由此又可证明：当采用不同的单元网格尺寸分析同一个问题，如果要使计算结果一致，那么参数取值就要与网格相关。
以上两组对比充分验证了经典裂缝带理论的网格相关性，而这种相关性可以进一步从构件的破坏模式进行解释。对于一根素混凝土梁，当跨中梁底出现第一根裂缝后，承载力就会急剧下降，其他位置就不会再出现新的裂缝，如图2(d)所示三个模型计算得到的变形和开裂应变云图(跨中挠度达到2 mm时)就很好地模拟了这种破坏特征。而进一步仔细观察各单元的应力-应变变化过程，当跨中截面单元率先达到混凝土轴心抗拉强度ft后开裂进入软化段本构时，其余截面单元则随着构件承载力的下降进行弹性卸载而无法进入开裂，因此无论单元如何划分，只有跨中截面的一列单元会进入开裂，一根裂缝总是在一个单元的范围内弥散，所以裂缝带宽hc自然应取为一个单元的尺寸lele。此外，图2(d)云图中的开裂应变数值也可印证这一规律，模型B和模型C的单元网格尺寸相比模型A缩小了一倍，而他们的开裂应变数值也相应放大了大约1倍，因为裂缝弥散成应变的标距缩小了，自然应变就放大了，可见裂缝弥散成应变的标距，也就是裂缝带宽hc，应该就是单元网格尺寸了。
值得关注的是，与上面这个算例类似，Bazant和 Oh[2]在论文中讨论的都是不配筋的素混凝土构件，用于模型验证的试验都是不配筋素混凝土的缺口试验，这些构件最大的特点是一旦在最薄弱的地方出现一条裂缝后，就不会在其他地方出现新的裂缝，自然具有显著的局部化特点。然而，实际工程结构中绝大部分都采用配筋混凝土构件，由于配筋的存在，混凝土的裂缝以某一个间距均匀地分布在构件中，这种裂缝的分布模式与素混凝土构件的集中分布模式截然不同。退一步讲，实际工程里即使采用素混凝土或配筋很少的混凝土构件，一旦达到混凝土轴心抗拉强度开裂后，就马上发生脆性破坏丧失承载力(如图 2(c)所示)，因此这种构件在实际工程中根本不允许出现裂缝，而此时再去研究如何计算开裂后的裂缝宽度则显得毫无意义。因此，针对素混凝土结构构件的经典裂缝带理论离结构工程实践具有一定的距离，只有将裂缝带理论根据配筋混凝土构件的力学特点进行改造和拓展，才能真正解决结构工程的问题，这也是接下来本文要进一步讨论的内容。
1.3 裂缝带理论在配筋混凝土构件中的表达形式
在配筋混凝土构件中，由于钢筋和混凝土之间的粘结作用，混凝土的裂缝以某一特定的间距分布出现，因此，在一个有限元网格 lele范围内有可能出现多根裂缝，这是与素混凝土构件在开裂行为上最显著的区别。基于这样的认识，可以将一个单元内的开裂应变看作这一个单元内所有裂缝宽度总和的弥散化，因此开裂应变εcr与一根裂缝的宽度w之间满足以下关系式：
式中，n为一个单元内的平均裂缝数量。
Bazant和Oh[2]定义开裂应变能gf如下：
图3所示为开裂应变和开裂应变能的定义。裂缝带模型认为，混凝土开裂后，其弹性拉应变εte随应力的降低而逐渐降低，其在开裂前积累的弹性应变能在开裂后是缓慢释放的，在受拉软化曲线上任意一点的应变ε，均可分解为弹性应变εte和开裂应变εcr，如图3(a)所示，式(5)所定义的开裂应变能gf表示的就是图3(a)中曲线下部涂色部分的面积。这里需要补充的是，一些通用有限元程序中的弥散裂缝模型与这一定义略有区别，譬如，MSC.MARC [11]认为混凝土开裂后，其弹性应变能瞬间释放，在受拉软化曲线上任意一点的应变ε就是开裂应变，则开裂应变能gf表示的面积比Bazant和Oh[2]定义的面积大，如图3(b)所示，这两个观点虽然对混凝土开裂机理的解释是不同的，但由于实际工程重点关注正常使用极限状态下(荷载水平大约为极限荷载的50%)裂缝稳定形成后的裂缝宽度值，此时这一差别就显得无关紧要，因为在混凝土从受拉到完全开裂(应力降为 0)的过程中发生的总开裂应变能 gf均指的是受拉应力-应变关系曲线的下部面积(如图1(b)所示)，且完全开裂点的开裂应变εcr均与总应变ε 相等。因此，裂缝开裂全过程的总开裂应变能 gf可写为：
式中，εtu为完全开裂点的总应变，如图3所示。
根据断裂能的定义，一根裂缝开裂全过程的断裂能Gf为：
式中，w0的含义为完全开裂(拉应力为 0)时的裂缝宽度，如图1(a)所示。
图3 开裂应变和开裂应变能的定义
Fig.3 Definitions of crack strain and crack strain energy
根据式(4)、式(6)和式(7)，可以推出裂缝开裂全过程gf与Gf的关系为：
根据平均裂缝间距lm的含义，即可得到一个单元内的平均裂缝数量n如下：
将式(9)代入式(8)，gf与Gf之间的关系可进一步写为式(10)，可以发现式(10)与单元网格尺寸lele无关。
由于 gf代表的为受拉应力-应变曲线的下部面积，若同样采用如图1(b)所示最简单的受拉线性软化形式时，gf可具体写为：
将式(11)代入式(10)，即可得到软化模量 Ets的表达式如下：
从式(12)可以得到一个非常重要的结论，那就是对于配筋混凝土构件，基于裂缝带理论推导得到的混凝土单轴应力-应变关系与单元网格尺寸无关，与素混凝土的结论截然相反。这一结论意味着当采用有限元方法进行配筋混凝土构件的裂缝宽度分析时，不必考虑单元网格划分对计算结果的影响，这给配筋混凝土构件的裂缝宽度分析带来了极大的便利，因为在实际结构的分析中，网格划分的选取往往受制于构件布置、计算代价等其他诸多因素，如果计算结果与网格无关，那么网格划分在具体操作过程中就会比较自由。
式(12)还表明，当混凝土材料一定时，断裂能Gf、混凝土开裂应力ft、混凝土弹性模量Ec都是确定的，则平均裂缝间距lm是影响混凝土受拉本构关系的最关键因素。裂缝越密，则软化模量 Ets就越小，曲线的下降段就越缓。这一规律同样可以从裂缝带理论的基本假定加以解释，裂缝越密集意味着单位长度内的裂缝数量越多，而单个裂缝开裂所需的断裂能Gf是一定的，越多的裂缝数量意味着构件受拉过程中耗散的能量越多，因此，构件因开裂导致的损伤退化发展也就越迟缓，表现在应力-应变关系上就是曲线下降段越缓。
最后，当完成非线性有限元分析后，可根据得到的开裂应变εcr按式(13)算得裂缝宽度值。同样的，这一公式也与网格尺寸无关。
分别对比式(12)与式(3)、式(13)与式(2)，对于配筋混凝土结构，裂缝带理论中的核心参数——裂缝带宽hc应取为平均裂缝间距lm，而非单元特征尺寸 lele，这从裂缝带的基本概念入手也非常容易解释，裂缝带表示一根集中的裂缝宽度均匀弥散成应变的范围，对于按一定间距lm分布的一系列裂缝，每一根裂缝恰好向两侧各 0.5lm范围内弥散才能恰好使所有裂缝弥散到开裂单元的所有范围内。
以上讨论均在理论层面，但在实践层面，混凝土开裂分析作为强非线性分析，常常会遇到难以回避的收敛性问题。软化模量 Ets对程序的收敛性具有十分关键的影响，许多时候，按照式(12)计算出来的软化模量Ets值较大(尤其是配筋率不太高的情况下，平均裂缝间距较大)，此时虽然理论上很完善，但输入到程序中却难以收敛，无法得到计算结果，那么这种理论上的完善就失去了意义。应对此问题时最常见的做法就是人为地调小软化模量 Ets值，使程序收敛，本文将该方法简称为“直接调整法”，如图1(c)所示，但这样的做法改变了材料的断裂能，从根本上违背了“断裂能是表征裂缝力学行为的本质参数”这一基本结论，其可能导致的误差将在后面的算例中做进一步的讨论。而另一种可行的做法是仍然坚持“采用断裂能来描述一根裂缝开裂全过程的力学行为”这一最基本的假定，在保证应力-应变关系曲线下包面积不变(也就是开裂应变能 gf不变)的情况下，通过人为调小开裂应力ft至ft,eq，来间接地调小软化模量，使程序便于收敛，本文将该方法简称为“间接调整法”，如图1(c)所示，这一做法的效果也将在后面的算例中进行讨论。
1.4 计算流程总结
在上述讨论的基础上，本文总结了采用通用有限元程序的弥散裂缝模型和分层壳单元计算裂缝宽度的流程，如图4所示。在这一计算流程中，平均裂缝间距lm占据了至关重要的位置，它架起了连接“应变”与“裂缝宽度”的桥梁，从而实现了有限元计算的“开裂应变”转换为“裂缝宽度”这一最终目标。关于该计算流程的具体说明如下：
步骤1：根据混凝土棱柱体抗压强度fc和最大骨料粒径Dmax，经过相应的强度换算采用式(1)计算断裂能Gf；
步骤 2：采用理论或经验公式计算平均裂缝间距lm；
步骤3：根据式(10)通过lm将断裂能Gf转换为开裂应变能gf，也就是应力-应变曲线的下包面积；
步骤4：由应力-应变曲线的下包面积gf，根据混凝土弹性模量Ec和混凝土轴心抗拉强度ft，采用式(12)算得受拉软化模量Ets，输入有限元程序进行非线性有限元分析；
步骤5：若程序不收敛，则将ft折减为ft,eq，重新按步骤4计算，若程序收敛，则提取有限元程序算得的开裂应变εcr；
步骤6：根据式(13)通过lm将开裂应变εcr转换为裂缝宽度w。
在有限元的理论体系里只存在“应变”的概念，但要完整地描述开裂行为必须要有“平均裂缝间距”和“裂缝宽度”两个物理量，因此，单纯依靠有限元方法，裂缝是不可解的，这就是有限元方法在处理非连续介质问题时的固有缺陷。这也从根本上解释了为什么在上述计算流程中，平均裂缝间距 lm需要游离于有限元裂缝分析过程之外预先单独计算确定，而无法在有限元裂缝分析过程之中自动解得。同时，这也表明平均裂缝间距在开裂应变与裂缝宽度之间的重要作用，确定平均裂缝间距的取值成为裂缝计算不可或缺的条件。
图4 采用通用有限元程序计算裂缝宽度的建议流程Fig.4 Proposed procedure for calculating the crack width using general-purpose finite element packages
目前，确定平均裂缝间距lm主要依赖于半经验半理论公式，通过试验数据回归是确定平均裂缝间距的最常用方法，此种方法在平均裂缝间距求解的机理上存在一定缺陷，更通用的计算方法仍有待进一步研究。
2 算例讨论
2.1 算例参数和模拟结果
为全面验证并讨论上述计算理论和流程，选取一由聂建国等[19]完成的承受负弯矩作用的简支组合梁静力加载试验进行有限元分析，分析的重点是定量预测混凝土板裂缝宽度的全过程发展。
图5所示为算例的基本信息及网格的划分。算例的具体参数和加载模式如图5(a)所示。基于通用有限元程序 MSC.MARC[11]中的弥散裂缝模型，采用梁-壳混合模型进行分析，钢筋混凝土板采用分层壳单元进行模拟(混凝土分为 5层)，钢梁采用纤维梁单元进行模拟(钢梁上翼缘、腹板、下翼缘的纤维划分数分别为16、20、16)。该梁-壳混合模型的单元类型、材料本构等建模信息详见文献[20—22]。分析采用位移收敛准则，误差限设为0.01。
平均裂缝间距lm是计算中的关键参数，由于承受负弯矩的组合梁其混凝土板近似处于轴心受拉状态，因此这里采用文献[23]建议的轴心受拉钢筋混凝土构件的平均裂缝间距lm计算公式，具体如下：
式中：c为保护层厚度；dr为钢筋直径；ρr为配筋率。
图5 承受负弯矩的简支组合梁[19]
Fig.5 Simply supported composite beam under hogging moment[19]
2.2 结果的网格依赖性
网格无关性是裂缝带理论推广到配筋混凝土构件时的重要特性，为验证这一结论，这里分别采用如图5(b)所示的4种不同的单元划分进行模拟，从模型A到模型D，单元网格由密到疏。
网格尺寸对模拟结果的影响如图 6所示。图6(a)给出了采用不同网格尺寸计算得到的跨中弯矩-跨中挠度曲线及其与试验结果的对比。不同单元网格的模型给出几乎完全一致的结果，说明结果确实与网格无关。所有的计算结果都与试验结果吻合良好，也验证了模型的准确性。此外，文献[20—22]对本算例采用的梁-壳混合模型开展了大量的试验验证，充分验证了该模型的可靠性。图 6(b)给出了跨中弯矩-跨中板分层壳单元顶层裂缝宽度曲线，网格尺寸对计算结果也几乎没有影响。此外，当混凝土板顶纤维达到开裂应力 ft后，计算曲线的裂缝宽度值发生突变，这也很好地印证了前述MSC.MARC[11]采用的弹性应变能集中释放的假定。
图6 网格尺寸对模拟结果的影响
Fig.6 Influence of mesh size on the simulation results
为进一步验证模型的网格无关性，图7对比了两个荷载水平下不同网格尺寸模型给出的混凝土板分层壳单元顶层开裂区范围，同样可以清楚地看到，不同的网格划分计算出的混凝土板开裂区范围也非常接近。以上从各个角度充分证明了单元网格的无关性，后续的分析均采用模型D的网格划分，进一步研究软化模量的调整方法和平均裂缝间距对裂缝宽度计算结果的影响。
图7 不同网格尺寸混凝土板开裂区范围计算结果对比Fig.7 Comparison of calculation results of cracking zone range of RC slab among different mesh sizes
2.3 软化模量相关收敛困难的应对策略
为应对基于裂缝带理论计算得到的软化模量Ets过大而导致的收敛困难，前述从理论层面讨论了两种应对策略：一种是直接人为调小软化模量 Ets值，即“直接调整法”；另一种是基于断裂能等效通过调小开裂应力ft来间接调小软化模量Ets值，即“间接调整法”。以下通过对本算例进行参数分析，进一步讨论这两种调整法的实际计算效果，以期对理论部分讨论的结论进行验证。
首先讨论“直接调整法”。软化模量Ets调整为裂缝带理论计算值的0.75倍和0.5倍，同时保持其他所有参数不变。采用“直接调整法”处理的计算结果如图8所示。图8(a)给出了跨中板分层壳单元顶层裂缝宽度随荷载的发展曲线对比，可见三根曲线有显著的差异。调小软化模量对开裂点以及刚开裂之后的一小段裂缝发展的计算结果没有影响，但随着荷载水平不断提高，裂缝不断发展，裂缝宽度的预测误差逐步增大，直到荷载水平大约超过极限承载力的50%后，裂缝宽度的预测误差才逐渐减小直至所有曲线又回到相同的结果。因此，“直接调整法”会导致正常使用荷载水平下显著的裂缝宽度预测误差，尽管如此，从图8(b)可以看到这一调整对宏观曲线的影响较小。
图8 “直接调整法”应对软化模量过大导致收敛困难的效果
Fig.8 Effect of direct adjustment method for overcoming the convergence difficulty due to large softening modulus
下面进一步讨论“间接调整法”，也是本文推荐的方法。因为软化模量的大小直接决定收敛的困难程度，为了使两种方法具有可比性，这里在保证开裂应变能不变的情况下通过调整开裂应力ft同样使软化模量Ets变为裂缝带理论计算值的0.75倍和0.5倍，即近似达到同样的收敛难易度。由式(11)可知，调整后的等效开裂应力 ft,eq与真实开裂应力ft的比值γ 按式(15)计算：
式中，β为调整后的软化模量与裂缝带理论计算得到的软化模量的比值，对本算例，当β 取为0.75和0.5时，γ 分别为0.878和0.726。
采用“间接调整法”处理的计算结果如图 9所示。图9(a)给出了采用“间接调整法”后跨中板分层壳单元顶层裂缝宽度随荷载的发展曲线对比，可见其误差情况与“直接调整法”完全相反，由于改变了开裂应力，开裂荷载以及刚开裂后的一小段裂缝发展有一定的误差，但很快随着荷载水平逐步提高(大约荷载水平提高到极限荷载的 20%之后)，三个模型几乎给出相同的裂缝宽度预测结果。因此，“间接调整法”在达到与“直接调整法”同样的软化模量折减效果的同时，也很好地保证了正常使用荷载水平下裂缝宽度的预测精度，同时由图9(b)还可以看到，“间接调整法”除了略微改变了开裂点附近的宏观曲线外，对整体宏观曲线的影响也很小，可见“间接调整法”是一种非常理想的应对软化模量过大导致收敛困难的策略。
图9 “间接调整法”应对软化模量过大导致收敛困难的效果Fig.9 Effect of indirect adjustment method for overcoming the convergence difficulty due to large softening modulus
为了更直观地展现两种调整法的计算效果，图10(a)和图 10(b)分别给出了两种调整法在不同荷载水平(定义为实际荷载与极限荷载的比值)下裂缝宽度值的预测误差，图10中w0为采用裂缝带理论得到的裂缝宽度值，是衡量误差的基准。从图 10中可以清楚地看到，“直接调整法”产生误差的主要荷载水平范围为20%~50%，“间接调整法”产生误差的主要荷载水平范围为20%以下，而实际混凝土结构构件正常使用阶段裂缝宽度验算最关心的荷载水平范围恰好落在 20%~50%，因此对于面向实际工程结构应用的裂缝宽度预测，应选用“间接调整法”。
图10 不同荷载水平下裂缝宽度的预测误差
Fig.10 Crack width prediction errors under different load levels
实际上，无论采用哪种调整法，在材料本构层面产生的误差主要都集中在材料完全开裂点(即拉应变达到图1(c)中的εtu)之前，从图8(a)和图9(a)中也可清楚地看到，当混凝土材料的拉应变达到完全开裂点εtu之后，裂缝宽度的预测误差就开始显著减小，而主要的误差都集中在材料达到完全开裂点εtu之前，因此要减小调整软化模量造成的误差，就需要尽可能提前εtu，使得材料完全开裂点尽可能避开混凝土结构构件正常使用状态的荷载水平，而图1(c)清楚地表明，要达到同样的软化模量折减程度，“间接调整法”求得的完全开裂点εtu相比“直接调整法”求得的完全开裂点有显著的提前，这也从基本原理上进一步解释了为什么应当选用“间接调整法”。
2.4 平均裂缝间距的影响
如前所述，平均裂缝间距lm是整个计算流程中的关键参数，它架起了连接“应变”与“裂缝宽度”的桥梁，因此作为讨论的最后一部分，本节将重点研究平均裂缝间距对计算结果的影响。将平均裂缝间距分别变化为式(14)计算值的0.5倍、0.75倍、1.25倍和1.5倍，其他参数均保持不变，代入图4所示的流程进行计算。图11(a)对比了跨中板分层壳单元顶层裂缝宽度随荷载的发展曲线，清晰地展示了平均裂缝间距对裂缝宽度计算结果的决定性作用，采用不同的平均裂缝间距得到的裂缝宽度结果差异很大，因此要得到准确可靠的裂缝宽度值，选取准确可靠的平均裂缝间距至关重要。而图11(b)表明，平均裂缝间距的改变对宏观曲线的影响仍然较小，可以忽略。
图11 平均裂缝间距对模拟结果的影响
Fig.11 Influence of the average crack spacing on the simulation results
3 关于网格无关性的其他佐证
“网格无关性”是本文研究的重要结论，为了进一步说明该结论的可靠性，以下给出了一些其他佐证材料供读者参考。
1) 代尔夫特科技大学(Delft University of Technology)的 Rots在他的博士论文[24]中对这个问题开展了一系列算例研究，发现当直接模拟钢筋与混凝土的粘结滑移效应时，裂缝计算结果和网格相关，而如果使钢筋与混凝土直接变形协调，裂缝宽度则与平均裂缝间距相关，他还引用了当年国际上很有代表性的 Collins等[25]完成的钢筋混凝土平板试验来证明这个结论。
2) de Borst曾于2002年在《Engineering Fracture Mechanics》中发表综述论文[26]，也发现了对于钢筋混凝土构件，在计算裂缝带宽时，需要考虑平均裂缝间距的影响。
3) 这里再补充一个算例来说明网格无关性结论，该算例为一单纯的钢筋混凝土简支板，如图12所示。计算结果如图 13所示，通过对比不同网格尺寸对应的跨中弯矩-跨中挠度曲线以及跨中弯矩-跨中板分层壳单元底层裂缝宽度曲线，同样可以得到网格无关性的结论。
图12 钢筋混凝土简支板
Fig.12 Simply supported reinforced concrete slab
4) 作为分层壳单元向一维化的退化，纤维梁单元近年来也广泛应用于结构分析，这种建模方式和分层壳单元一样不直接模拟钢筋与混凝土层之间的粘结滑移，笔者团队通过算例研究也证明了采用纤维梁模型模拟钢筋混凝土构件裂缝宽度同样与网格无关，详见文献[27]。
图13 网格尺寸对模拟结果的影响
Fig.13 Influence of mesh size on the simulation results
4 结论
本文讨论了采用通用有限元程序的弥散裂缝模型和分层壳单元计算钢筋混凝土裂缝宽度的理论基础与方法。本文首先简要回顾了基于弥散裂缝模型的经典裂缝带理论，并结合某无筋素混凝土梁的有限元算例讨论了该理论的局限性。在此基础上，从配筋混凝土结构构件的裂缝发展特点出发，推导了裂缝带理论在配筋混凝土构件中的表达形式，讨论了数值收敛问题的应对策略，并给出了计算流程。最后，以某承受负弯矩的简支组合梁混凝土板开裂分析为例，全面验证了本文理论讨论部分的关键结论。总结本文研究，可得如下重要结论：
(1) 对于无筋素混凝土结构构件，单裂缝软化的局部化效应导致有限元计算结果与网格密切相关，裂缝带理论中的裂缝带宽应取为单元特征尺寸。
(2) 对于配筋混凝土结构构件，多根裂缝按照平均裂缝间距均匀分布，有限元计算结果与网格无关，裂缝带理论中的裂缝带宽应取为平均裂缝间距。
(3) 平均裂缝间距将有限元分析中的“应变”概念和工程设计中的“裂缝宽度”概念紧密联系起来，是本文建议计算流程中最重要的参数，对计算结果具有决定性的影响。
(4) 当遇到软化模量过大而导致的数值收敛问题时，建议在保持断裂能等效的前提下同时调整混凝土轴心抗拉强度和软化模量值。
在本文的讨论中，“平均裂缝间距”始终是最核心的关键词，因此，更通用化的平均裂缝间距计算方法非常值得进一步研究，作为初步的探索，笔者所在团队尝试将断裂能判据引入有限元方法求解平均裂缝间距，相关讨论详见文献[28]。此外，“受拉刚化效应”本文未讨论，因为要考虑该效应不可避免需对通用有限元程序进行二次开发，故不在本文的讨论范围内，笔者在文献[12]和文献[27]中对此问题有专门的讨论，读者可查阅。
[1] GB 50010-2010, 混凝土结构设计规范[S]. 北京: 中国建筑工业出版社, 2010.GB 50010-2010, Code for design of concrete structures[S]. Beijing: China Architecture & Building Press, 2010.(in Chinese)
[2] Bazant Z, Oh B. Crack band theory for fracture of concrete [J]. Material and Structures, 1983, 16(3): 155-177.
[3] 陶慕轩, 樊健生, 聂建国, 等. 型钢混凝土柱-钢桁梁组合节点抗震性能理论分析[J]. 工程力学, 2009,26(11): 152-160, 171.Tao Muxuan, Fan Jiansheng, Nie Jianguo, et al,Theoretical analysis on aseismic behavior of steel reinforced concrete column-steel truss beam composite joints[J]. Engineering Mechanics, 2009, 26(11): 152-160, 171. (in Chinese)
[4] 聂建国, 王宇航. ABAQUS中混凝土本构模型用于模拟结构静力行为的比较研究[J]. 工程力学, 2013, 30(4):59-67, 82.Nie Jianguo, Wang Yuhang. Comparison study of constitutive model of concrete in ABAQUS for static analysis of structures [J]. Engineering Mechanics, 2013,30(4): 59-67, 82. (in Chinese)
[5] 何树凯, 王来贵. 钢筋混凝土结构开裂分析在 ANSYS中的实现[J]. 辽宁工程技术大学学报, 2007(增刊 2):125-127.He Shukai, Wang Laigui. Realization of RC structure dehiscence analysis with ANSYS [J]. Journal of Liaoning Technical University, 2007(Suppl 2): 125-127. (in Chinese)
[6] 陈力, 方秦, 还毅, 等. 对 ABAQUS中混凝土弥散开裂模型的静力特性分析[J]. 解放军理工大学学报(自然科学版), 2007(5): 478-485.Chen Li, Fang Qin, Huan Yi, et al. Analysis on static performances of smeared cracking model for concrete in ABAQUS [J]. Journal of PLA University of Science and Technology, 2007(5): 478-485. (in Chinese)
[7] 张娟霞, 唐春安, 周秀艳, 等. 基于高性能计算的钢筋混凝土构件等间距裂缝形成过程研究[J]. 工程力学,2009, 26(3): 161-167.Zhang Juanxia, Tang Chun’an, Zhou Xiuyan, et al. Study on fracture spacing formation process of reinforced concrete specimen based on high performance calculation[J]. Engineering Mechanics, 2009, 26(3): 161-167. (in Chinese)
[8] 任重翠. 钢筋混凝土剪力墙拉剪性能试验研究[D]. 北京: 中国建筑科学研究院, 2018.Ren Chongcui, Experimental study on tension-shear performance of reinforced concrete shear wall [D].Beijing: China Academy of Building Research, 2018. (in Chinese)
[9] Namdar A, Darvishi E, Feng X, et al. Effect of flexural crack on plain concrete beam failure mechanism: A numerical simulation [J]. Frattura ed Integrità Strutturale,2016, 10(36): 168-181.
[10] Lang L, Zhu Z M, Zhang X S, et al. Investigation of crack dynamic parameters and crack arresting technique in concrete under impacts [J]. Construction and Building Materials, 2019, 199: 321-324.
[11] MSC. MARC 2012 [M]. Santa Ana, CA: MSC Software Corp, 2012.
[12] Xu L Y, Nie X, Zhou M, Tao M X. Whole-process crack width prediction of reinforced concrete structures considering bonding deterioration [J]. Engineering Structures, 2017, 142: 240-254.
[13] Wittmann F H, Slowik V, Alvaredo A M. Probabilistic aspects of fracture energy of concrete [J]. Materials and Structures, 1994, 27(9): 499-504.
[14] Cornelissen H A W, Hordijk D A, Reinhardt H W.Experimental determination of crack softening characteristics of normalweight and lightweight concrete[J]. Heron, 1986; 31(2): 45-56.
[15] 过镇海, 张秀琴. 砼受拉应力-变形全曲线的试验研究[J]. 建筑结构学报, 1988, 9(4): 45-53.Guo Zhenhai, Zhang Xiuqin. Experimental investigation of complete stress-deformation curves of concrete in tension [J]. Journal of Building Structures, 1988, 9(4):45-53. (in Chinese)
[16] van Mier J G M, Schlangen E, Vervuurt A. Tensile cracking in concrete and sandstone: Part2 – Effect of boundary rotations [J]. Materials and Structures, 1996,29(3): 87-96.
[17] Hurbut B. Experimental and computational investigation of strain-softening in concrete [D]. Boulder, Colo;University of Colorado, 1985.
[18] CEB-Comite Euro-International du Beton-EB-FIP Model Code 1990 [S]. Bulletin D Information No.213/214,Lausanne, 1993.
[19] 聂建国, 李法雄, 樊健生, 等. 不同翼板形式组合梁受弯性能试验研究[J]. 公路交通科技, 2009, 26(10): 53-58.Nie Jianguo, Li Faxiong, Fan Jiansheng, et al.Experimental study on flexural behavior of composite beams with different styles of concrete flange[J]. Journal of Highway and Transportation Research and Development, 2009, 26(10): 53-58. (in Chinese)
[20] 聂建国, 陶慕轩. 多高层钢-混凝土组合框架结构体系弹塑性分析模型[J]. 建筑结构学报, 2010, 31(7): 1-12.Nie Jianguo, Tao Muxian. Model for elasto-plastic analysis of multistory and highrise steel-concrete composite frame systems [J]. Journal of Building Structures, 2010, 31(7): 1-12. (in Chinese)
[21] Nie J G, Tao M X, Cai C S, et al. Modeling and investigation of elasto-plastic behavior of steel–concrete composite frame systems [J]. Journal of Constructional Steel Research, 2011, 67(12): 1973-1984.
[22] 陶慕轩. 钢-混凝土组合框架结构体系的楼板空间组合效应[D]. 北京: 清华大学, 2012.Tao Muxian. Slab spatial composite effect of steel-concrete composite frame structural systems [D].Beijing: Tsinghua University, 2012. (in Chinese)
[23] 叶列平. 混凝土结构[M]. 第2版, 上册. 北京: 清华大学出版社, 2005.Ye Lieping. Concrete structures [M]. 2nd ed, Volume Ⅰ.Beijing: Tsinghua University Press, 2005. (in Chinese)
[24] Rots J G. Computational modeling of concrete fracture[D]. Delft: Delft University of Technology, 1988.
[25] Bhide S B, Collins M P. Reinforced concrete elements in shear and tension [D]. Publication 87-02, University of Toronto, Department of Civil Engineering, 1987.
[26] de Borst R. Fracture in quasi-brittle materials: a review of continuum damage-based approaches [J]. Engineering Fracture Mechanics, 2002, 69(2): 95-112.
[27] Xu L Y, Nie X, Tao M X. Rational modeling for cracking behavior of RC slabs in composite beams subjected to a hogging moment [J]. Construction and Building Materials, 2018, 192: 357-365.
[28] Wang J J, Tao M X, Nie X. Fracture energy-based model for average crack spacing of reinforced concrete considering size effect and concrete strength variation [J].Construction and Building Materials, 2017, 148: 398-410.