基于離散宏單元的砌體結構高效非線性分析方法

郭 勇，余丁浩，李 鋼

(大連理工大學海岸和近海工程國家重點實驗室，遼寧，大連 116024)

砌體結構因其造價低廉、施工方便、性能優異而廣泛存在于工業與民用建筑中，此類結構的主要材料為塊體和砂漿組成的非均質復合材料，其力學行為由塊體、砂漿和二者之間粘結界面的力學性能共同確定。塊體和砂漿為典型的脆性材料，抗拉強度遠小于抗壓強度。二者之間粘結界面的力學性能差，界面受拉行為與材料粘性屬性有關，界面剪切行為同時與材料粘性和摩擦屬性有關。在極端環境荷載作用下，砌體結構將表現出顯著的非線性特征。

目前，砌體結構非線性分析模型按照建模方式的不同可分為分離式和整體式兩種。分離式模型將砂漿、塊體及二者粘結界面分開建模，各部分的材料性質通過描述各自力學性能的本構關系表征，模型精細化程度高，但計算代價大。Louren?o和Rots[1]在塊體-砂漿粘結界面模型中采用復合屈服準則考慮砌體墻的拉、壓、剪非線性行為，并在非線性行為中考慮了材料的塑性軟化特征；Sandoval 等[2]利用該模型完成了多孔黏土磚砌體墻和灌漿配筋砌體墻[3]的精細化數值模擬，與試驗結果的對比分析結果顯示該模型能夠反映砌體的真實破壞形態； Louren?o[4]和Roca 等[5]對砌體結構非線性分析的數值計算問題進行了綜述，認為分離式模型最突出的問題在于計算代價高昂，嚴重地限制了其在大型工程結構中的應用；Minga等[6]針對分離式模型應用于大型砌體結構計算代價高昂的缺點，采取將結構分解為若干部分，分配給不同處理器并行計算的方式進行非線性分析。整體式分析模型將塊體和砂漿假定為連續的勻質體，建模過程簡單，非線性分析所需計算量小，此類模型屬于宏觀分析模型，可以較好地平衡結構分析的計算效率和精度，其關鍵在于建立能夠合理反映結構實際震害特點和關鍵失效模式的簡化建模策略，當前已有眾多學者開展了此方面研究，并提出了多種整體式分析模型，如勻質化模型[7-9]、等效框架模型[10-18]等，Grande 等[19]在梁單元的端部設置剪切界面，提出了可考慮砌體墻肢剪切失效的梁單元模型，將所提模型與節點區剛性單元聯合使用模擬砌體墻的彎曲和剪切非線性行為；Addessi 等以Timoshenko 梁理論為基礎，分別在分布式非線性梁單元模型[20]和集中非線性梁單元模型[21]中設置塑性剪切鉸，并用該塑性鉸模擬墻肢的剪切失效機制；Caliò等[22]提出了一種考慮砌體墻面內主要失效模式的平面離散宏單元模型，該模型用非線性彈簧構造剪切單元和無厚界面單元，通過耦合兩種單元的變形模擬砌體墻的彎曲和剪切非線性行為，目前，該模型經不斷改進被應用于無筋砌體墻在靜力作用下[23]和動力作用下[24]的平面外受力特性分析、填充墻鋼筋混凝土框架結構的平面內抗震分析[25-26]和平面外抗震分析[27-28]、無筋砌體結構的地震易損性分析[29]等，研究結果均表明該模型具有良好的分析精度和計算效率。綜上所述，精細化的分離式模型適用于結構構件的細部分析和模擬，而整體式模型在大型砌體結構的分析和模擬方面具有更大優勢。

無論整體模型或分離模型，通常采用傳統變剛度法進行非線性求解，不可避免地需對切線剛度矩陣進行實時更新和分解，這導致隨問題規模增大計算效率急劇降低。為了避免結構大規模整體切線剛度矩陣的實時更新和分解，國內外學者利用土木工程結構材料非線性的局部化特征，發展出一系列高效數值分析方法，例如多尺度方法、子結構方法、重分析方法等。其中，多尺度有限元法通過構造多尺度基函數，將局部非線性區域的精細化模型與其余部位的宏觀模型進行有效耦合，避免了傳統有限元法采用單一小尺度建模，大幅縮減結構自由度，Drosopoulos 和Stavroulakis[30]、 Giambanco 等[31]、 Silva 等[32]和Krejci 等[33]采用多尺度方法對砌體結構進行了非線性分析，分析結果證明該方法計算效率明顯優于傳統單一尺度建模方法。Clough 和Wilson[34]提出子結構方法，該方法依據邊界力和變形協調條件將整體結構劃分為線彈性子結構和非線性子結構，并通過靜力凝聚的方式消去線彈性子結構的內部自由度，大幅降低結構自由度數目，提高了計算效率，但該方法需要預判非線性子結構的劃分區域，且彈、塑性子結構一旦劃分就無法改變，因而適用性有限。為克服其局限性，孫寶印和古泉等[35-36]對整體結構的單元彈塑性狀態進行實時判斷，僅將進入非線性狀態的構件作為隔離塑性子結構，使大規模的非線性分析問題轉化為整體結構的線彈性分析和規模較小的局部隔離子結構的非線性分析。重分析方法是一類利用初始結構計算信息快速求解修改后結構響應的計算方法，能夠避免結構局部修改后進行耗時的完全分析，此類方法最初用于結構設計，隨后被眾多學者用來提升局部材料非線性分析的計算效率。重分析方法主要包含直接法和近似法，直接重分析法[37]將局部結構修改或局部非線性變形產生的剛度矩陣表示為初始彈性剛度矩陣的低秩修正形式，并基于Shearman-Morrison-Woodbury 公式進行問題求解，避免了整體切線剛度矩陣耗時的實時更新和分解，組合近似法[38]是基于預處理共軛梯度法發展而來的一種高效的近似重分析方法，該方法通過對變化的剛度矩陣進行降階處理來縮減問題的計算規模，可快速求解大范圍非線性結構的響應。擬力法[39-41]將構件的變形分為彈性和塑性兩部分，可在控制方程中保持初始彈性剛度矩陣不變并通過僅更新小規模塑性剛度矩陣實現局部非線性問題的高效求解，受其概念啟發，Li等[42-44]提出了隔離非線性法，該方法通過對非線性材料的應變進行分解，并與有限元理論結合，可將切線剛度矩陣轉化為初始彈性剛度矩陣的低秩修正形式，結合Woodbury 公式實現了一般局部非線性有限元分析問題的高效求解，為克服該方法隨非線性區域規模擴大時效率降低的不足，Li等[45-46]通過引入組合近似法的思想對該方法使用的Woodbury 求解公式進行了深度優化。上述學者針對結構局部非線性特征發展了多種高效數值分析方法，極大推動了結構分析領域的發展。

地震作用下砌體結構的非線性變形狀態往往集中于樓梯間、縱橫墻交接處等局部區域，具有典型的局部化特征，因此，利用其局部非線性特征發展適用的高效分析方法對提升計算效率具有積極意義。本文基于整體式空間離散宏單元模型，并充分利用結構局部非線性特征，提出了一種砌體結構高效非線性分析方法，首先建立剪切單元和無厚界面單元，并通過耦合兩者之間的變形模擬砌體失效模式，隨后，結合隔離非線性思想將剪切單元等效斜向彈簧的軸向變形和無厚界面單元上下表面的相對變形分解為線彈性和非線性兩部分，并通過附加的塑性自由度描述其中的非線性部分，從而在控制方程中實現了線性與非線性剛度矩陣的分離，以該控制方程為基礎將結構切線剛度矩陣表示為初始彈性剛度的低秩修正形式，最后，通過引入Woodbury 公式求解控制方程，使非線性分析過程的主要計算量集中于對一個小規模的非線性矩陣進行更新與分解，有效提升了非線性分析效率。數值模擬與試驗結果的對比驗證了本文方法的準確性，基于時間復雜度理論的計算效率分析結果表明該文方法在非線性分析效率方面較傳統變剛度法具有明顯優勢。

1 離散宏單元模型

Caliò等[22]基于整體式建模策略提出了砌體結構的空間離散宏單元模型，如圖1 所示，其根據砌體墻窗洞構造將墻體離散成一系列剪切單元和無厚界面單元，剪切單元為由四個剛性板和兩個對角斜向非線性彈簧構成的鉸接四邊形結構，無厚界面單元由一系列不同功能的非線性彈簧構成。不同剪切單元之間或剪切單元與支座之間通過無厚界面單元相互連接，砌體墻的平面內、外失效模式通過耦合剪切單元和無厚界面單元的變形進行模擬。

圖1 空間離散宏單元模型Fig. 1 Spatial discrete macro-element model

砌體墻斜截面剪切破壞是受損砌體結構中最普遍的失效模式，當墻體在豎向荷載和水平地震作用下產生的主拉應力超過墻體的抗拉強度時，平行于地震作用方向的墻體產生斜向裂縫，如圖2(a)所示。該破壞模式通過剪切單元進行模擬，剪切單元將墻肢等效為發生純剪變形的均質板，非線性行為通過兩個對角斜向非線性彈簧進行描述，如圖3(a)所示。彈簧的初始彈性剛度可基于墻肢與均質板線彈性變形范圍內的能量等效[22]確定，屈服強度可根據Turnsek-Cacovic 準則[22]確定，如圖4(a)所示為斜向彈簧的單軸滯回模型[22,24]，模型中考慮了材料屈服軟化和卸載剛度退化等非線性行為，卸載剛度退化系數可取為β=0.8[24]。

圖2 砌體墻平面內外主要失效模式Fig. 2 Main in-plane and out-of-plane failure mechanisms of a masonry portion

圖3 離散宏單元模型模擬砌體墻平面內外主要失效模式(為簡化說明，圖(b)、(c)、(d)僅畫出一根斜向彈簧)Fig. 3 Simulation of the main in-plane and out-of-plane failure mechanisms of a masonry portion by means of the macro-element (For simplicity, only one diagonal spring is shown in figures (b), (c) and (d))

砌體墻正截面彎曲破壞模式表現為墻體受拉區的開裂和受壓區的壓碎，如圖2(b)所示。該破壞模式通過在墻肢截面處劃分的若干勻質彈塑性纖維模擬，纖維的受力狀態通過界面單元中m行n列法向非線性彈簧進行描述(圖3(b))，在多個剪切單元與界面單元組合而成的砌體墻模型中，這些界面彈簧僅影響與其相連的剪切單元的一半區域(如圖1 所示，墻體左上豎向界面單元僅影響單元左右兩側陰影區域)，圖4(b)為這些彈簧的單軸滯回模型[22,24]，模型中考慮了材料卸載剛度退化行為，卸載剛度退化系數取為 β=0.8[24]。

當砌體墻豎向荷載較小或砌體與砂漿界面摩擦系數較低時，墻體易發生正截面剪切滑移破壞，表現為沿砌體的灰縫出現貫通截面的滑動摩擦行為，如圖2(c)所示。為此，沿界面單元縱向設置一根剪切非線性彈簧描述墻體的正截面剪切滑移破壞，如圖3(c)所示，該剪切彈簧的本構關系采用理想剛塑性模型[23-24](本文數值計算時將圖4(c)所示理想彈塑性模型的彈性剛度取為大數)，屈服強度根據Mohr-Coulomb 準則[22]確定，模型中未考慮材料的卸載剛度退化行為，即 β=0。

圖4 砌體材料的單軸本構關系Fig. 4 Uniaxial constitutive laws of masonry material

當砌體結構縱、橫墻間連接薄弱時，在沿橫墻方向的水平地震作用下結構內外墻交接處易產生豎向裂縫，此時縱墻向外甩出發生平面外彎曲、剪切和扭轉變形，此時墻體變形如圖2(d)所示。模型中沿界面單元橫向設置兩根剪切非線性彈簧(圖3(d))模擬墻體的平面外剪切、扭轉變形，每根彈簧的單軸滯回模型如圖4(c)所示，屈服強度根據Coulomb 失效準則確定[23-24]。

上述剪切單元通過無厚界面單元連接建立砌體墻片模型，如圖1 所示，進一步將墻片模型通過共結點的方式連接，并根據結構樓蓋形式及支承條件設定邊界約束，即可建立整體結構分析模型，進而實現整體結構破壞形式和損傷演化過程的模擬。

2 非線性隔離宏單元模型

基于隔離非線性法基本思想，將材料應變或單元變形分解為線彈性和非線性兩部分，圖5 以單軸材料的非線性本構關系為例給出了應變分解的基本原理，假設在某個增量步內，材料的非線性狀態由B點到C點，根據初始彈性模量Ee將材料應變增量分解為兩部分：

圖5 非線性材料應變分解Fig. 5 Decomposition of nonlinear material strain

式中： Δε′為材料的線彈性應變增量； Δε′′為材料的非線性應變增量。與應變增量 Δε相對應的應力增量 Δσ可表示為：

式(2)表明，無論何種本構模型，應變分解后任意時刻的材料應力均可表示為恒常的初始彈性模量與線彈性應變的乘積。

2.1 剪切單元模型

選擇剪切單元的形心作為基本結點，每個基本結點k有7 個自由度，其中6 個為剛體運動自由度，1 個為剪切變形自由度，如圖6 所示，其中oxyz為單元局部坐標系。在迭代過程的任一線性化增量步中，基本結點k位移增量向量可表示為：

如圖6 所示，單元每個剛性板的中心設置1 個外部結點，共4 個外部結點。其中：①號、②號結點沿局部x軸布置；③號、④號結點沿局部z軸布置。每個外部結點r有6 個自由度，即有6 個廣義位移和6 個廣義力，其中僅在xz面內的剪力使單元發生剪切變形。外部結點r位移增量向量為：

圖6 剪切單元結點自由度Fig. 6 Shear element nodal degrees of freedom

圖7 給出了單元的剪切變形模式(考慮小變形假設)，將③、④號外部結點的連線定義為單元剪切變形模式的主線，該主線在xz面內的轉動變形Δγy(基本結點剪切變形自由度)將引起主線端部外部結點③、④的平動位移dz·Δγy和其他外部結點①、②的轉動位移 Δγy。

圖7 剪切單元剪切變形模式Fig. 7 Shear deformation mode of shear element

在單元剪切變形模式下，等效斜向彈簧的變形如圖8 所示，根據圖中幾何關系，可建立剪應變 Δ γy和彈簧變形 Δd之間的關系：

圖8 等效斜向彈簧變形圖Fig. 8 Deformation of equivalent diagonal spring

式中： Δd為剪切單元斜向彈簧變形增量； Φdiag為變形矩陣。

基于圖5 所示的隔離非線性方法應變分解思想，可將等效斜向彈簧變形 Δd分解為線彈性和非線性兩部分：

式中，線彈性部分 Δd′等于單元等效斜向彈簧初始彈性剛度的倒數乘以等效斜向彈簧內力增量ΔN，即：

非線性部分 Δd′′等于彈簧變形增量 Δd與線彈性變形增量 Δd′的差值。斜向彈簧內力增量可表示為其初始彈性剛度與線彈性變形增量的乘積，即：

式中，Ddiag,e為剪切單元的等效斜向彈簧的初始彈性剛度。

2.2 無厚界面單元模型

基于多點線性約束方法(multi-point constraints,MPC)，剪切單元每一外部結點r的位移可以用基本結點k剛體運動引起的位移和剪切變形(圖7 所示剪切變形模式)引起的轉動位移或平動位移表示。在局部坐標系中，求解外部結點r的位移時，定義以基本結點k為起點、外部結點r為終點的向量=(dx,dy,dz)(如圖9 給出了k=2，r=④2時的向量，④2表示基本結點2 的④號外部結點)，則單元外部結點r的位移函數為：

圖9 無厚纖維界面單元Fig. 9 Zero-thickness fiber interface element

式中，等號右端第一項為基本結點k剛體運動引起的位移項，第二項和第三項分別為剪切變形引起的轉動位移項和平動位移項，其中 Δγrot和Δγdisp分別為剪切變形模式下控制外部結點r轉動位移和平動位移的轉動自由度，erot=(,,)和edisp=(,,) 分 別為 Δγrot和 Δγdisp的單位方向向量，在圖7 所示的單元剪切變形模式下，控制外部結點①、②位移的轉動自由度 Δγrot=Δγy，Δγdisp=0，erot=(0,1,0)，控制外部結點③、④位移的轉動自由度 Δγrot=0 ， Δγdisp=Δγy，edisp=(0,1,0)。

圖9 表示出無厚纖維界面單元連接基本結點1 的③號外部結點(即圖中③1結點)和基本結點2 的④號外部結點(即圖中④2結點)。從圖9 可知相鄰③1、④2號外部結點的相對位移產生了界面相對變形，基于此可將界面單元相對位移函數表示為：

式中： Δuo={ΔuΔvΔwΔθxΔθyΔθz}T為界面單元相對變形增量；I6為6 階單位矩陣； ΦI為變形矩陣。

將無厚界面單元的相對變形 Δuo分解為線彈性和非線性兩部分：

為方便后文對圖3(b)、圖3(c)、圖3(d)描述的單元變形分量展開討論，對界面單元相對變形增量進行初等變換后寫為：

式中： ΔuCB={ΔwΔθyΔθx}T為單元法向彈簧模擬壓彎受力行為的變量； ΔuinS={Δu}T為局部x向切向彈簧控制面內剪切滑移變形的變量；ΔuoutST={ΔvΔθz}T為局部y向切向彈簧模擬面外剪扭受力行為的變量。基于此可將分解后的界面單元變形寫為如下形式：

式(16)中，DCB需通過對單元中各法向彈簧的彈性剛度集成得到，即：

式中：ECB和A分別為彈簧的彈性模量和影響面積；xi和yi為第i根彈簧的局部坐標，彈簧沿局部y向布置m行，沿局部x向布置n列，m×n為彈簧的數量，DinS根據單元中沿局部x向切向彈簧的彈性剛度得到：

式中：EinS代表彈簧的彈性模量；沿單元局部x向設置一根切向彈簧，影響面積AinS為單元的面積；DoutST需對單元中沿局部y向兩根彈簧的彈性剛度集成得到：

式中：EoutST和A分別為彈簧的彈性模量和影響面積；xj表示第j根彈簧的局部坐標。

1) 根據隔離非線性法中的變形分解原理，單元中第i根法向彈簧的相對位移增量可分解為：

實際分析中采用Newton-Raphson 迭代法求解，理想彈塑性模型會導致剛度矩陣奇異而出現求解失敗，為了改善數值計算的收斂性，本文將本構關系中屈服后剛度與彈性剛度的比值取為0.001。當彈簧本構關系采用理想彈塑性模型時，根據圖5 所示非線性區段內的材料應變分解思想，相應的彈簧應力增量可表示為：

對單元中各法向彈簧的應力增量進行集成，可得與 ΔmCB對應的單元法向內力和彎矩表達式：

式中： ΔNz為單元法向內力； ΔMy和 ΔMx為單元彎矩。將式(22)代入式(16)中可建立與各法向彈簧的相對位移的關系：

對比式(23)和式(12)可知，單元相對變形分量 ΔuCB和非線性相對變形分量可分別表示為：

式(25)表明，當單元中所有法向彈簧的非線性相對位移均為零值時，單元的相對非線性變形分量才為零向量(此時相應塑性自由度處于未激活狀態)，若部分彈簧的非線性相對位移為非零值，則無論這些彈簧數量的多寡，單元非線性相對變形分量均為非零向量(此時相應的塑性自由度被激活)。

2) 沿單元局部x向僅設置一根切向彈簧，與ΔminS對應的單元剪力為：

式中， ΔVx為單元剪力。根據單元線彈性變形表達式(12)可知當單元中該切向彈簧處于彈性狀態時，單元非線性相對變形分量為零向量(相應塑性自由度未激活)，若為非零值，則為非零向量(相應塑性自由度被激活)。

3) 同樣地，對單元中沿局部y向切向彈簧的剪切相對位移進行分解，將兩根彈簧的應力增量進行集成，得與 ΔmoutST對應的單元剪力和扭矩為：[ΔVyΔMz]T=

式中： ΔVy為單元剪力； ΔMz為單元扭矩。進一步根據單元線彈性變形表達式(12)可知單元相對變形分量 ΔuoutST和非線性相對變形分量分別為：

式(29)表明，當單元中兩根切向彈簧的非線性剪切相對位移均為零值時，單元的非線性相對變形分量才為零向量(相應塑性自由度未激活)，若其中一根彈簧的非線性剪切相對位移為非零值，單元非線性相對變形分量為非零向量(相應塑性自由度被激活)。

2.3 單元控制方程

上述剪切單元和界面單元的變形分解后，單元內力增量與其變形增量和非線性變形增量相關，式(5)和式(9)、式(10)分別給出了剪切單元和界面單元變形增量的計算模型，單元非線性變形通過 Δd′′和進行描述，其中每個元素代表一個塑性自由度。

圖9 建立的兩基本結點單元可作為模型的基礎單元，包含兩個剪切子單元及一個界面子單元，建立其虛功方程如下：

式中： δu為單元的虛變形向量；ΔT=[ΔNT]T為單元內力； δa為單元虛結點位移向量； Δf為結點荷載增量。將單元位移函數式(5)、式(9) 、式(10)和內力增量表達式(8)、式(16)代入虛功方程式(30)，整理可得：

此外，考慮單元中各彈簧的非線性本構關系可建立補充方程，對于進入非線性狀態的單元內力使用下式求解：

將式(32)代入式(31)，可得基礎單元的隔離非線性控制方程：

式中：ke為單元的初始彈性剛度矩陣；k′為描述單元結點力與非線性變形之間聯系的系數矩陣；為描述單元非線性行為的系數矩陣。具體通過下式計算：

式中，B矩陣為代表基本結點自由度與兩剪切子單元斜向彈簧變形 Δd和界面子單元相對變形 Δuo之間關系的變形矩陣。對B矩陣，可先將式(9)代入式(10)得到界面單元的變形矩陣，再將式(5)擴充為與其統一階數的矩陣形式，最后進行合并得到。De=diag(DI,e,Ddiag,e)，且Ddiag,e=diag(,)，和分別為圖9 中相鄰1、2 號剪切單元的等效斜向彈簧的初始彈性剛度，Dt為相應的切向剛度矩陣。

3 結構控制方程求解

將基礎單元控制方程進行集成，可得整體結構的隔離非線性控制方程：

式中： ΔX和 ΔF分別為整體結構的結點位移增量向量和結點荷載增量向量；為整體結構的非線性變形增量向量，其中僅包含了結構中呈激活狀態的m個塑性自由度；Ke為整體結構的初始彈性剛度矩陣，其規模為n×n階，階數等于結構位移自由度數n；K′為描述整體結構結點力與非線性變形之間聯系的系數矩陣，其規模為n×m階；為描述整體結構材料非線性行為的系數矩陣，其規模為m×m階。當結構發生局部材料非線性變形時，激活的塑性自由度數目m將遠小于結構自由度數目n(即m＜＜n)。

式(38)表明：當m遠小于n時，任意時刻結構的切線剛度矩陣可表示為初始彈性剛度矩陣的低秩攝動形式。

采用完全Newton-Raphson 迭代法進行結構非線性分析，引入Woodbury 公式(式(39))在每個迭代步中對相應式(38)進行高效求解：

式中，KSchur為Schur 補矩陣：

從式(39)可以看出，結構整體剛度項Ke在整體分析過程中始終保持初始彈性狀態，僅需在分析前進行一次分解即可，當材料非線性集中在結構的局部區域時，控制方程中塑性自由度數目遠小于結構自由度數目(即m＜＜n)，此時Woodbury公式在求解結構切向位移反應時的主要計算開銷將僅集中于小規模的Schur 補矩陣(m×m階)的更新分解[43]，從而可避免傳統變剛度法中對整體切線剛度矩陣(n×n階)的反復合成與分解運算，實現計算效率的大幅提升。

4 數值算例

4.1 試驗驗證

為了驗證本文提出的砌體結構非線性分析方法的正確性，以Anthoine 等[47]的擬靜力試驗數據為依據，對2 片不同高度的磚砌體墻進行數值模擬。墻體的材料屬性和幾何屬性分別如表1 和表2所示，2 片墻體在豎向壓應力下的低周反復荷載試驗結果如圖10 所示，其中圖10(a)為Wall-1 正截面彎曲破壞的試驗滯回曲線，可以看出，曲線捏縮較為顯著，墻體耗能能力較弱但強度在變形較大時沒有明顯退化；圖10(b)為Wall-2 斜截面剪切破壞的試驗滯回曲線，曲線的包圍面積較大，墻體耗能能力相對較強，但存在嚴重的剛度和強度退化，試驗詳細信息見參考文獻[47]。

表1 墻體材料屬性Table 1 Material properties of masonry walls

表2 墻體幾何屬性Table 2 Geometry properties of masonry walls

圖10 砌體墻片試驗滯回曲線Fig. 10 Experimental behaviour of simple piers

采用本文方法模擬2 片砌體墻在豎向荷載和水平荷載共同作用下的面內滯回行為，每個模型使用1 個剪切單元和1 個設置于底端的無厚界面單元。模型中各類非線性彈簧的力學參數取值根據表1 中砌體墻力學參數確定。兩個模型數值模擬結果如圖11 所示，分別為Wall-1 和Wall-2 基底剪力與頂點位移之間的滯回曲線，此外，圖12給出了Caliò和Pantò等求解的Wall-1 和Wall-2 兩片砌體墻宏單元模型在低周反復荷載作用下的滯回曲線(CP 方法求解結果)[22]。對比砌體墻數值模擬與試驗結果：三者剛度和強度整體上吻合良好，其中CP 方法求解Wall-1 滯回曲線捏縮更為顯著，這是由模型中非線性彈簧本構關系中滯回規則的差異所導致的。為更充分驗證本文方法的正確性，圖13 提取了三組滯回曲線的骨架曲線，從圖中知，對比Wall-1 頂端最大位移處荷載P1，本文方法求解結果和CP 方法求解結果與試驗結果的相對誤差分別為12.38%和16.14%；對比Wall-2 頂端最大位移處荷載P2，本文方法求解結果和CP 方法求解結果與試驗結果的相對誤差分別為9.51%和25.73%。可見：本文方法求解誤差小于CP 方法，本文所提方法可以有效地模擬砌體墻的失效模式，驗證了該分析方法的正確性。

圖11 砌體墻片本文方法計算滯回曲線Fig. 11 Numerical simulation results of piers(proposed method)

圖12 砌體墻片CP 方法計算滯回曲線Fig. 12 Numerical simulation results of piers (CP method)

圖13 砌體墻片骨架曲線Fig. 13 Skeleton curves of simple piers

4.2 動力非線性分析

為進一步對本文方法進行驗證，對某4 層約束砌體結構進行地震反應分析，該結構平面簡圖如圖14 所示。結構層高為3 m，總高度為12 m，縱墻上窗洞的尺寸為1.5 m ×1.8 m，門洞的尺寸為1.0 m ×1.8 m。在結構縱橫墻交接處設置構造柱(如圖14 所示)，結構所有縱橫墻在樓層標高處設置鋼筋混凝土圈梁，構造柱和圈梁采用C20 混凝土，截面尺寸為250 mm ×250 mm，縱筋采用4 根HRB335直徑16 mm 的鋼筋；墻體厚度為250 mm，材料為燒結普通磚，強度等級為MU10，混合砂漿強度等級為M10。結構樓面恒載為5.49 kN/m2，活載為2.0 kN/m2，屋面恒載為6.371 kN/m2。設防烈度為8 度(0.2g)，Ⅱ類場地，地震動選取汶川地震中臥龍波(幅值為400 cm/s2)，加速度時程如圖15 所示。圖16 為本文方法建立的該結構數值分析模型，其中，模型假設在結構首層底端設置固定約束，圈梁和構造柱采用隔離非線性纖維梁單元進行模擬，共劃分1894 個隔離非線性纖維梁單元[48]，砌體墻劃分為4198 個剪切單元和9770 個無厚界面單元，該空間模型共有39 934 個位移自由度。圈梁與構造柱和砌體墻的相互連接通過在接觸面設置界面單元的方式實現，另外，在樓板的每一區格設置兩根斜向剛性桿傳遞水平荷載和協調樓面結點位移。為簡化模擬，本例中僅考慮纖維界面單元的軸向和彎曲非線性變形，其剪切和扭轉變形假設為彈性。

圖14 結構平面簡圖Fig. 14 Plane of the masonry structure

圖15 臥龍地震波加速度記錄Fig. 15 Acceleration record of Wolong earthquake

圖16 約束砌體結構數值分析模型Fig. 16 Discrete macro-element model of the masonry structure

結構動力響應分析之前對結構進行模態分析，發現結構第一陣型是一階Y向平動，一階陣型的周期為T1=0.14 s，第二陣型是一階X向平動，二階陣型的周期是T2=0.1 s。分析結構在Y向水平地震作用下的動力響應，結構外縱墻與外橫墻交接處各樓層(圖16 示)位移時程響應如圖17所示，結構的基底剪力時程響應如圖18 所示，根據各樓層位移時程曲線發現該結構底層層間位移角最大，峰值達到0.6%，依據文獻[49]中約束砌體結構0.4%的倒塌限值判定該砌體結構已達到倒塌的破壞狀態。本砌體結構主要震害集中在底層，外橫墻墻肢的損傷分布如圖19 所示，墻體下部1 層、2 層的剪切單元均發生剪切破壞，下部損傷明顯比上部嚴重。

圖17 各樓層位移時程曲線Fig. 17 Time-history curves of all floors

圖18 結構基底剪力時程曲線Fig. 18 Base shear force curve of the masonry structure

圖19 外橫墻損傷分布圖Fig. 19 Damage distribution of outer transversal wall

圖20 給出了塑性自由度數隨增量步的變化曲線，在3.265 s(增量步為653 步)之前結構處于彈性狀態，塑性自由度數為0。非線性分析過程中產生的最大塑性自由度數為6545，占總自由度數的16.4%，而每個增量步的平均塑性自由度為1997，僅占總自由度數的5%，說明結構具有顯著的局部非線性特征，采用式(39)進行求解時，由于需進行實時更新和分解的Schur 補矩陣階數等于塑性自由度數目，局部非線性引發的較少的塑性自由度可保證本文方法實現大幅度的效率提升。

圖20 塑性自由度時程曲線Fig. 20 Time-history curves of plastic degrees of freedom

為了進一步論證本文方法的高效性，將傳統非線性分析方法(矩陣分解采用LDLT法)和本文方法的時間復雜度[50]進行統計對比，如圖21 所示，傳統非線性分析方法和本文方法的平均時間復雜度分別為1.54 ×1010和1.40 ×109，傳統方法的平均時間復雜度是本文方法的11 倍，證明了本文方法的高效性。

圖21 時間復雜度時程曲線Fig. 21 Time-history curves of time complexity

5 結論

本文通過耦合剪切單元和無厚界面單元來模擬砌體結構失效行為，以隔離非線性法理論為基礎，通過將剪切單元中等效斜向彈簧軸向變形和無厚界面單元上下表面的相對變形分解為線彈性和非線性兩部分，推導出砌體結構離散宏單元的隔離非線性控制方程，采用Woodbury 公式求解控制方程，得到如下結論：

(1) 本文利用砌體結構材料非線性的局部化特征建立了結構非線性隔離宏單元模型，該模型基于多點線性約束方法使每個基本結點的自由度大幅縮減，相比于已有宏觀模型和精細模型，其在大型砌體結構非線性分析時既能保證較高的計算精度也能使模型規模維持在相對較低的程度。同時，模型中激活的塑性自由度數目等于結構控制方程中小規模非線性矩陣的維度，其可直觀地衡量砌體結構中產生增量非線性變形的區域規模。

(2) 通過分離結構彈塑性狀態使砌體結構非線性分析時僅對小規模非線性矩陣進行更新和分解，可以避免傳統變剛度法中對整體切向剛度矩陣的更新與分解，該方法的時間復雜度明顯低于傳統變剛度法，可以實現砌體結構非線性問題的高效求解。