基于極限分析的三葉墻抗壓強度預測模型研究

蔣宇洪，楊 娜

(北京交通大學土木工程系，北京 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)所得抗壓強度預測值的準確性顯著高于常見的三葉墻抗壓強度預測公式，對古建筑保護具有積極意義。