基于极限分析的三叶墙抗压强度预测模型研究

蒋宇洪，杨 娜

(北京交通大学土木工程系，北京 100044)

砌体结构是历史最为长远、应用最为广泛的复合材料结构之一，随着砌筑材料、砌筑工艺的变换，会呈现不同的力学性能。过去砌筑高大的砌体结构时，为了满足承载能力和稳定性的需求，工匠选择大厚度墙体作为主要承重构件。受限于当时的材料加工、运输手段，很难采用巨型砌块砌筑厚墙，于是工匠提出了一种利用多种规格砌块的砌筑工艺：如图1 所示，采用较大且规整的砌块砌筑墙体的两片外层，外层所夹空间填充着夯土或较小石块-砂浆的混合物。这样的结构被称为三叶墙(three-leaf wall)，常出现在城墙、教堂、宫殿、民居等古代砌体建筑中；墙体外层被称为外叶墙(outer-leaf)，填充内层被称为内叶墙(innerleaf)。

图1 三叶墙(西藏古建石砌体)Fig.1 Three-leaf wall (Tibetan ancient stone masonry)

目前，国外学者们对三叶墙进行了许多试验研究，着重于探索构造、材料、加固手段对三叶墙力学性能和行为的影响。Binda 等[1]对不同外-内叶墙连接形式的三叶墙进行抗压试验，指出具有键状连接的三叶墙仅比平整连接的三叶墙的强度有小幅的提高。Silva 等[2]、Corradi 等[3]、Valluzzi等[4]和Oliveira 等[5]分别对未加固和已加固的三叶墙进行抗压试验，对比分析了不同加固手段对三叶墙承载能力的影响。

砌体结构数值模拟的收敛性一直是个难题，而三叶墙构造比普通砖墙构造更为复杂，给数值模拟提出了更大的挑战。为了规避有限元法的收敛难题，同时达到更快的计算效率，一些学者采用极限分析对砌体结构的极限强度进行研究，极限分析常用于土体稳定问题的求解[6-8]，但在复合材料结构中的应用较少。Sloan[9]提出了一种三角形单元的有限元塑性极限分析下限法，提供了对结构划分任意网格并进行极限分析的基础。近年来，一些外国学者将其用于砌体结构的极限强度预测。Sutcliffe 等[10]提出将Sloan[9]的有限元塑性极限分析下限法运用于砌体结构剪力墙的抗剪强度计算；Milani 等[11]基于极限分析法，假定砌块为刚体，砂浆传递的应力为多项式分布，提出了一个分析砌体结构面内力学性能的模型，并研究了砌体构造参数和材料性质对结果的影响。Casolo等[12]在Milani 等[11]所提模型的基础上，在厚度方向上进行了拓展，对三叶墙的面外抗弯性能进行了研究。

目前利用极限分析对砌体结构进行的研究很少，而且几乎没有关于三叶墙承载能力的研究。本文基于砌体结构有限元塑性极限分析下限法，将受压三叶墙破坏模式与材料试验、小型砌体试件试验的结果相结合，赋予了本构模型参数。在此基础上，考虑了受压三叶墙的不均匀压力荷载分布，提出了一个三叶墙抗压强度预测模型。根据文献的试验数据，利用该模型预测了多个三叶墙的抗压强度，并与试验结果和抗压强度预测公式的结果进行了对比，说明了该模型的准确性。

1 有限元塑性极限分析下限法

如果一结构内部的应力场均满足平衡微分方程，不违背屈服条件，同时边界应力场满足边界条件，这样的应力场称为静力许可应力场[13]。下限极限分析法就是在所有静力许可应力场对应的荷载中，找到最大的荷载，即为极限荷载。该方法的基本假设包括：1)材料简化为理想刚塑性；2)结构濒临破坏前的变形仍为小变形；3)材料满足Drucker 公设。

本文以Sutcliffe 等[10]的砌体结构有限元塑性极限分析下限法为基础，将结构离散成三角形单元和界面单元，把单元的节点应力作为未知量，根据平衡条件、应力连续条件、边界条件、屈服准则得到约束方程，由外荷载和节点应力的关系式获得目标函数，最后基于约束方程和目标函数，提出线性规划问题，求得外荷载的极限值。值得注意的是，有限元塑性极限分析下限法与常用的有限元法有很大差别，该方法基本变量是节点应力，不涉及单元应变和结构变形，所以没有要求满足几何方程、变形协调方程和物理方程。此外，该方法所施加的边界条件固定不变，因此只能对单调静力加载问题进行求解。

1.1 单元离散和平衡条件

有限元塑性极限分析下限法将结构离散成如图2 所示的三角形单元，三角形单元之间利用零厚度的应力间断界面单元进行连接。图2 中的x向代表水平方向；y向代表竖向；界面单元与x轴正向的夹角为θ。

图2 有限元离散Fig.2 Finite element discretization

三角形单元上的每个节点上均有三个平面应力分量，同时设定单元内部的应力分量是线性分布的。于是该三角形单元的任意一点的应力可用节点应力进行表示，如式(1)所示：

式中：{σ}是三角形上任意一点的平面应力，{σ}={σxσy τxy}T；N1、N2、N3为三角形单元的形函数，其表达式见文献[14]；{σj}为三角形节点j的节点应力，{σj}={σx jσyj τxy j}T。

任意一点的应力需要满足平衡微分方程，即需满足式(2)、式(3)。

式中，γ为材料重度。

将式(1)代入式(2)、式(3)并整理，即可得到平衡条件对应的约束方程Ae，见式(4)：

1.2 应力连续条件

如图2 所示，单元边界上的应力可用应力分量 σx、σy、τxy表示，也可用法向应力 σn和切应力τ表示，两组应力分量可用式(5)、式(6)转换：

有限元塑性极限分析下限法不同于常见的有限法，相邻单元的节点可在同一坐标重合，利用应力间断界面单元相连。为了满足应力连续，允许重合节点的 σx、σy、τxy不等，但必须保证重合节点的 σn和τ相等。对于如图2 所示的两单元，当重合节点①、节点②的边界应力满足式(7)时，即代表满足应力连续条件。

式中：σn1、σn2分别是节点①、节点②的法向应力；τ1、τ2分别是节点①、节点②的切应力。

由于单元应力是线性分布的，当作为线段端点的节点应力保持连续时，该线段上任意一点的应力也保持连续。将式(5)、式(6)代入式(7)并整理，得到应力连续条件对应的约束方程Ai，见式(8)。

1.3 边界条件

实际问题中，结构边界往往受到法向、剪切向的分布荷载，如图3 所示。分布荷载和其合力的关系式，如式(9)所示：

图3 应力边界条件Fig.3 Stress boundary condition

式中：q、t分别为边界上法向和切向分布荷载；P、Q分别是法向和切向分布荷载合力；S为结构的边界。

因为单元应力是线性分布的，式(9)可写为式(10)：

如图3 所示，为了简化模型，假设处于同一条边界上的节点①、节点②的应力满足：σn1=σn2，τ1=τ2；将式(5)、式(6)代入式(10)中，经过整理，可以得到约束方程Ab，见式(11)。

值得注意的是，当外荷载为线性分布荷载或不均匀分布荷载时，可以通过细分边界单元，对边界施加等效的多段均布荷载代替。

1.4 屈服准则

本文基于Lourenco 等[15]提出的改进摩尔-库仑屈服模型，对摩尔-库仑模型加入了最大拉应力和压应力的限制，如图4 所示。为了满足线性规划的使用限制，对该屈服面进行线性近似简化。同时，将该屈服面分割成三部分，分别表示不同的破坏形式：帽盖面、剪切破坏面、立断面分别反映的是材料受压破坏、剪切破坏、受拉破坏。

图4 摩尔-库仑模型屈服面Fig.4 Yield surface of Mohr-Coulomb model

为了引出该模型的屈服方程，先给出经典的摩尔-库仑屈服方程，如式(12)所示：

式中，c、φ分别为材料的内粘聚力和内摩擦角。

将式(5)、式(6)代入式(12)，即可得到利用平面应力分量表示的屈服方程，再代入不同大小的φk、ck，即可分段表示线性近似屈服函数，如式(13)所示：

式中，下标k代表不同段屈服面的编号，k为1、2、3 时，分别对应图4 屈服面上半部分的立断面、剪切破坏面、帽盖面，k为4、5、6 时，分别对应下半部分的立断面、剪切破坏面、帽盖面。

对式(13)进行整理，可用式(14)表示。

该屈服方程的具有多个不同大小的φk、ck。屈服面代表剪切破坏时，φk、ck与材料真实的内粘聚力内摩擦角φ、c相等；屈服面代表抗压破坏和抗拉破坏时，φk是根据所需屈服面形状设定的，参考Sutcliffe 等[10]的设定值，设定φ1=φ4=89°，φ3=φ6=150°，ck可结合抗压、抗拉强度得到。综上，φk、ck可按式(15)设定：

式中：fc、ft分别为材料的抗压、抗拉强度。

将式(14)进行整理，可得屈服准则的约束方程Ay，见式(16)：

由于单元应力线性分布，当单元的各个节点应力满足屈服条件时，则单元中任一点的应力均不会违背屈服条件。

1.5 目标函数

对于结构极限荷载问题，需要在保证满足约束条件的同时，找到能使外荷载尽量大的应力场。将式(5)和式(6)代入式(10)，可以得到外荷载和节点应力的关系式，如式(17)所示：

1.6 下限线性规划问题

通过综合上述函数约束函数和目标函数，可以获得如式(18)所示的线性规划问题。通过对该问题的求解，所得C{σj}l为结构极限荷载的下限解，此时的{σj}即为极限荷载对应的结构整体应力解。

式中，C=[c1c2]。

2 三叶墙抗压强度预测模型

本文在Sutcliffe 等[10]的砌体结构有限元塑性极限分析下限法的基础上，结合常见受压三叶墙破坏模式，基于材料试验和小型砌体试件试验赋予了本构模型参数。在此基础上，考虑了不均匀压应力边界条件，对三叶墙抗压强度问题进行求解。

2.1 受压三叶墙破坏模式

三叶墙的受压破坏通常呈现脆性，并且破坏时往往出现较为明显的破坏现象。本文基于文献[1 -5,16 -17]中的三叶墙抗压试验的破坏现象，总结了如图5 所示的常见受压三叶墙破坏模式：

图5 受压三叶墙破坏模式Fig.5 Failure mode of three-leaf wall under compression

1)灰缝受压而膨胀，导致砌块受拉破坏：砌块出现较为密集的竖向裂缝，裂缝可能贯穿多个砌块和灰缝。砌体受到垂直于灰缝轴线的压力作用时，由于砌块和灰缝的弹性模量、泊松比差异，灰缝往往比石头横向膨胀得更多。为保持砌块和灰缝之间的位移连续性，接触界面粘结力和摩擦力使得灰缝受到侧向位移约束。这导致灰缝处于三轴压应力状态，砌块同时处于单轴压缩和双向受拉的状态。最终，砌块的拉应力达到了砌块的抗拉强度，砌块上最终出现受拉裂缝。

2)边界约束和压力作用下，砌块受剪破坏：处于边界的砌块出现了由外向内的斜竖向裂缝。在靠近边界的部分，砌块同时受到竖向压力和水平向剪切力，砌块切应力超过砌块抗剪强度而破坏。

3)内叶墙压坏：内叶墙内部出现多条竖向裂缝。内叶墙是由比较碎小的砌块和砂浆组成，承载能力很小；当内叶墙的压应力达到了其抗压强度后，出现大量受压裂缝。

4)接触界面开裂：外-内叶墙的接触界面出现裂缝，可能是外、内叶墙分离，也可能是该区域的砂浆开裂。三叶墙在承受压力时，内叶墙和外叶墙的接触界面会发展横向拉应力，当拉应力超过界面抗拉强度时，界面开裂。

5)外叶墙受压失稳：内叶墙膨胀，界面开裂并出现明显的局部空隙，外叶墙砌块向外鼓出并有脱落的趋势。外叶墙在承受压力的同时，内叶墙的膨胀使得外叶墙产生了面外变形，该外叶墙受压失稳并丧失大部分承载能力。

从破坏现象可以看出：受压三叶墙破坏模式1)～破坏模式4)均是不同组成部分的材料破坏；而破坏模式5)可以视为一个有侧向位移的受压柱的结构破坏，该破坏与外叶墙外鼓变形、试件长细比等因素直接相关。而本文采用的极限分析法不涉及变形变量，因此，本文的三叶墙抗压强度预测模型仅能反映除外叶墙受压失稳以外的受压破坏模式1)～破坏模式4)。

上述破坏模式中，没有灰缝压应力超过灰缝材料抗压强度，导致丧失承载能力的破坏。这是因为多数灰缝材料在受到超过自身抗压强度的压力之后，往往被压缩密实[18]，如图6 所示。但该灰缝仍然能起到传递压应力的作用，不会导致砌体结构承载能力丧失。

图6 压缩密实的灰缝[18]Fig.6 Compressed joint[18]

2.2 本构模型参数

除了常见的砌块、灰缝砂浆的基本材性试验，学者们针对三叶墙特性，设计了特定的小型砌体试件试验。对于外叶墙，部分学者[1-2,17]模仿外叶墙构造，砌筑如图7(a)所示的外叶墙试件，将该试件的抗压强度作为外叶墙抗压强度，该方法可以较大程度地还原外叶墙的构造，但是外叶墙试件所用砌块往往是等厚度的，不能还原外叶墙砌块厚度有变化的情况。而Oliveira 等[5]将少数砌块和砂浆砌成如图7(b)所示的棱柱体试件，将该棱柱体的抗压强度作为外叶墙的强度。该方法的构造相似程度较低，但是用料较少，可多次试验，降低结果离散性。对于内叶墙，大多数学者[1,5,17]将碎石和砂浆混合，制成如图7(c)所示的圆柱体，把该圆柱体的抗压强度作为内叶墙强度。

图7 三叶墙小型砌体试件Fig.7 Small masonry specimen of three-leaf wall

砌体领域内，学者们往往对材料试件进行抗压试验和抗折试验，分别得到该材料的抗压强度和抗折强度。然后根据抗折强度与抗拉强度比例为1.5∶1 的关系，得到该材料的抗拉强度。而本文屈服准则所用的材料参数不仅需要材料的抗压、抗拉强度，还需要材料的内粘聚力和内摩擦角。本文借鉴Milani 等[11]的转换方法，抗压强度fc、抗拉强度ft与内粘聚力c、内摩擦角φ的关系式如式(19)、式(20)所示。

将式(19)和式(20)联立方程求解，得式(21)、式(22)，将所得的抗压强度和抗拉强度代入，即可获得内粘聚力c、内摩擦角φ：

对于三叶墙的不同组分，均可采用式(14)、式(15)作为屈服函数，代入不同的材料参数即可。对于砌块屈服模型，代入砌块材性试验所得参数，可以反映三叶墙破坏模式2)；对于内叶墙屈服模型，代入内叶墙圆柱体抗压试验所得参数，可以反映三叶墙破坏模式3)；对于灰缝的屈服模型，因砌体抗压破坏基本不受灰缝材料抗压强度的影响，采用Lourenco 等[15]的方法，代入外叶墙试件或棱柱体的抗压强度作为灰缝的抗压强度，其余参数采用灰缝材料的参数，可反映三叶墙破坏模式1)；对于外-内叶墙界面屈服模型，采用灰缝材料的参数，反映的是三叶墙破坏模式4)。

2.3 不均匀压应力分布的边界条件

对于三叶墙抗压试验，外叶墙和内叶墙构造本身存在显著的差异，导致竖向刚度不同，压应力分布不均匀。为了得到三叶墙合理的承载能力预测，本文借鉴Egermann 等[16]和Toumbakari 等[19]对三叶墙压应力的分析，提出一个考虑三叶墙不均匀压应力的边界约束条件。

顶部平均压应力分布如图8 所示，边界平均应力、应变定义式见式(23)。

图8 三叶墙顶部压应力分布Fig.8 Distribution of compressive stress on top of three-leaf wall

假设忽略了泊松比的影响，竖向平均压应力和平均压应变的关系式可表达为式(24)。

在砌体结构抗压试验中，往往通过大刚度的加载梁对试件施加压力。那么，外、内叶墙的竖向平均压应变是相等的，即，代入式(24)，推导可得式(25)：

将式(25)进行整理，得到三叶墙顶部压应力的约束方程At，见式(26)。At是考虑了不均匀应力分布的边界条件约束方程，因此，在规划问题中，At是Ab的一部分。

2.4 三叶墙有限元模型

本文模型采用了如图9 所示的简化分离式建模方法：对于外叶墙部分，灰缝利用界面单元代替，砌块利用扩大化的砌块单元代替，保证砌体结构整体尺寸不发生改变；对于内叶墙部分，简化为均匀统一的各向同性材料；对于外-内叶墙接触界面，采用界面单元代替。

图9 三叶墙简化分离式模型及网格Fig.9 Simplified micro model and mesh of three-leaf wall

如图9 所示，本文借鉴李泽[20]的单元划分方法，对三叶墙模型进行有限元网格的划分，将结构划分成一系列小方格，一个小方格由四个三角形网格组成。在三叶墙受压问题中，结构和边界条件沿竖向中线对称，因此可以将结构简化为如图10 所示的半结构：原本为对称轴的半结构边界设定为滑动支座，切应力为0；墙体侧面通常不受外力作用，设置该半结构左侧法向应力、切应力均为0；同时，墙体顶部受到不均匀的压应力作用，具体设置方法见本文2.3 节。

图10 三叶墙半结构边界条件及网格Fig.10 Boundary conditions and mesh of semi-structure of a three-leaf wall

本文采用MATLAB 编制了主体计算程序，利用软件自带的linprog 函数，对该线性规划问题进行求解。本文方法的流程图如图11 所示。

图11 三叶墙抗压强度预测模型流程图Fig.11 Flow chart of prediction model of compressive strength of three-leaf walls

3 结果分析

3.1 算例选择

目前对三叶墙的抗压性能已经进行了许多试验研究，但因为砌筑材料、砌筑工艺的多样性，导致不同类型三叶墙的结果差异较大。外叶墙砌块的尺寸规格是影响砌体结构不均匀性和性能的重要因素。如图12 所示，常见的三叶墙砌块种类包括：毛石、毛料石、细料石、黏土砖等，砌块形状规则程度、表面平整程度从前到后逐渐变高，对应试件的不均匀程度和试验结果离散程度逐渐变低[3]。

图12 三叶墙砌块种类Fig.12 Type of block of three-leaf walls

考虑文献所能提供材料参数的全面程度，本文选取了Binda 等[1]、Silva 等[2]、Oliveira 等[5]、Demir 等[17]的试验，进行抗压强度的计算。表1是文献三叶墙的材料性能参数，基本取自于原文献试验所得数据。

表1 三叶墙抗压试验的材料参数和结果Table 1 Material parameters and results of compressive tests of three-leaf walls

3.2 经典三叶墙抗压强度预测公式

为了说明本文所预测的三叶墙抗压强度的准确性，分别采用本文方法与抗压强度预测公式，对三叶墙抗压强度进行预测。以下是学者们[1,4,17]曾提出的一些三叶墙抗压强度的预测公式及相应假定：

1)外叶墙的刚度远大于内叶墙的刚度，因此，外叶墙承担了绝大部分的荷载。

2)假设外、内叶墙根据各自的体积承担荷载。

3)假设外、内叶墙根据各自的体积和调整系数承担荷载。

式中：fw为三叶墙抗压强度；fo、fi分别为外、内叶墙强度；θo、θi分别为外、内叶墙的抗压强度调整系数；Vo、Vi分别是外、内叶墙的体积。

值得注意的是，本文不进行式(29)的对比，因为该公式的准确性很大程度的依赖于调整系数θo、θi的大小。目前该系数设定往往取决于学者的经验判断，缺乏合理力学解释，导致预测结果不够客观。

3.3 三叶墙抗压强度预测结果的分析

为了证明本文抗压强度预测模型的准确性，接下来将对模型所得的抗压强度和应力场进行分析。图13 是分别利用本文方法和抗压强度预测公式所得的抗压强度。由计算结果可以看出：式(27)、式(28)往往高估墙体的抗压强度。这是因为预测公式假定了三叶墙的破坏模式是外叶墙压坏或外、内叶墙同时压坏，没有考虑到其他破坏模式的发生；此外，该公式没有考虑到三叶墙压应力分布不均的状态。对比本文方法所得强度与试验结果，除了Silva 等[2]试验的结果误差为37%，其他试验误差均在8%～20%范围内，说明了本文抗压强度模型的有效性。Silva 等[2]的毛石砌块三叶墙强度预测值的误差较大，是因为毛石砌块不规则的规格导致外叶墙试件与三叶墙中的真实外叶墙存在构造差距 。

图13 预测抗压强度与试验结果的对比Fig.13 Comparison of predicted and test values of compressive strength

图14(a)是本文模型所得的竖向应力云图，该云图是以半结构应力解为基础，考虑单元应力的线性分布特性，利用轴对称特性获得。图14(b)是采用分离式建模，对三叶墙进行弹性分析所得的应力云图[22]。可以看出图14(a)和图14(b)中的应力云图分布规律相近，这说明了本文模型可以获得三叶墙结构较为合理的应力场分布。压应力云图表明：三叶墙处于一个不均匀受压的状态，外、内叶墙协同承受竖向压力荷载时，刚度更小的内叶墙有下凹趋势，为了平衡相邻内叶墙区域的变形趋势，导致外叶墙内侧压应力更大，内叶墙外侧压应力更小。

图14 压应力分布对比分析/MPaFig.14 Comparison of compressive stress distributions

3.4 网格细化对预测结果的影响

为了探讨单元网格细化程度对三叶墙预测抗压强度的影响，本文以Demir 等[17]的试验试件为算例，建模方法与本文前述内容相同，指定了如图15所示的6 个由粗到细的网格划分形式，相应的三叶墙抗压强度计算结果如表2 所示。由表2 可以看出，本文方法的计算精度与网格的精细程度密切相关，网格设置越精细，计算精度越准确，但计算效率也会变得更慢。当单元网格精细到一定程度后，计算开始结果趋于稳定。

图15 Demir 等[17]试验的有限元网格划分Fig.15 Finite element meshes of test by Demir et al[17]

表2 不同网格划分所得抗压强度比较 /MPaTable 2 Comparison of compressive strength obtained with different meshes

4 结论

本文以有限元塑性极限分析下限法为基础，根据材料试验、小型砌体试件试验结果，并结合受压三叶墙破坏模式，赋予了本构模型参数。同时考虑压应力分布不均匀的边界条件，提出了三叶墙抗压强度预测模型。主要结论如下：

(1)该模型规避了常见有限元法的收敛问题，计算高效，可以基于小型试验结果，对三叶墙的抗压强度作出预测。

(2)所得抗压强度预测值的准确性显著高于常见的三叶墙抗压强度预测公式，对古建筑保护具有积极意义。