基于離散宏單元的砌體結(jié)構(gòu)高效非線性分析方法

郭 勇，余丁浩，李 鋼

(大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧，大連 116024)

砌體結(jié)構(gòu)因其造價(jià)低廉、施工方便、性能優(yōu)異而廣泛存在于工業(yè)與民用建筑中，此類結(jié)構(gòu)的主要材料為塊體和砂漿組成的非均質(zhì)復(fù)合材料，其力學(xué)行為由塊體、砂漿和二者之間粘結(jié)界面的力學(xué)性能共同確定。塊體和砂漿為典型的脆性材料，抗拉強(qiáng)度遠(yuǎn)小于抗壓強(qiáng)度。二者之間粘結(jié)界面的力學(xué)性能差，界面受拉行為與材料粘性屬性有關(guān)，界面剪切行為同時(shí)與材料粘性和摩擦屬性有關(guān)。在極端環(huán)境荷載作用下，砌體結(jié)構(gòu)將表現(xiàn)出顯著的非線性特征。

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

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

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

1 離散宏單元模型

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

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

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

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

圖3 離散宏單元模型模擬砌體墻平面內(nèi)外主要失效模式(為簡(jiǎn)化說明，圖(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))

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

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

圖4 砌體材料的單軸本構(gòu)關(guān)系Fig. 4 Uniaxial constitutive laws of masonry material

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

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

2 非線性隔離宏單元模型

基于隔離非線性法基本思想，將材料應(yīng)變或單元變形分解為線彈性和非線性兩部分，圖5 以單軸材料的非線性本構(gòu)關(guān)系為例給出了應(yīng)變分解的基本原理，假設(shè)在某個(gè)增量步內(nèi)，材料的非線性狀態(tài)由B點(diǎn)到C點(diǎn)，根據(jù)初始彈性模量Ee將材料應(yīng)變?cè)隽糠纸鉃閮刹糠郑?/p>

圖5 非線性材料應(yīng)變分解Fig. 5 Decomposition of nonlinear material strain

式中： Δε′為材料的線彈性應(yīng)變?cè)隽浚?Δε′′為材料的非線性應(yīng)變?cè)隽?。與應(yīng)變?cè)隽?Δε相對(duì)應(yīng)的應(yīng)力增量 Δσ可表示為：

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

2.1 剪切單元模型

選擇剪切單元的形心作為基本結(jié)點(diǎn)，每個(gè)基本結(jié)點(diǎn)k有7 個(gè)自由度，其中6 個(gè)為剛體運(yùn)動(dòng)自由度，1 個(gè)為剪切變形自由度，如圖6 所示，其中oxyz為單元局部坐標(biāo)系。在迭代過程的任一線性化增量步中，基本結(jié)點(diǎn)k位移增量向量可表示為：

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

圖6 剪切單元結(jié)點(diǎn)自由度Fig. 6 Shear element nodal degrees of freedom

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

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

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

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

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

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

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

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

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

2.2 無厚界面單元模型

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

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

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

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

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

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

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

式中： ΔuCB={ΔwΔθyΔθx}T為單元法向彈簧模擬壓彎受力行為的變量； ΔuinS={Δu}T為局部x向切向彈簧控制面內(nèi)剪切滑移變形的變量；ΔuoutST={ΔvΔθz}T為局部y向切向彈簧模擬面外剪扭受力行為的變量?；诖丝蓪⒎纸夂蟮慕缑鎲卧冃螌憺槿缦滦问剑?/p>

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

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

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

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

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

實(shí)際分析中采用Newton-Raphson 迭代法求解，理想彈塑性模型會(huì)導(dǎo)致剛度矩陣奇異而出現(xiàn)求解失敗，為了改善數(shù)值計(jì)算的收斂性，本文將本構(gòu)關(guān)系中屈服后剛度與彈性剛度的比值取為0.001。當(dāng)彈簧本構(gòu)關(guān)系采用理想彈塑性模型時(shí)，根據(jù)圖5 所示非線性區(qū)段內(nèi)的材料應(yīng)變分解思想，相應(yīng)的彈簧應(yīng)力增量可表示為：

對(duì)單元中各法向彈簧的應(yīng)力增量進(jìn)行集成，可得與 ΔmCB對(duì)應(yīng)的單元法向內(nèi)力和彎矩表達(dá)式：

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

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

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

2) 沿單元局部x向僅設(shè)置一根切向彈簧，與ΔminS對(duì)應(yīng)的單元剪力為：

式中， ΔVx為單元剪力。根據(jù)單元線彈性變形表達(dá)式(12)可知當(dāng)單元中該切向彈簧處于彈性狀態(tài)時(shí)，單元非線性相對(duì)變形分量為零向量(相應(yīng)塑性自由度未激活)，若為非零值，則為非零向量(相應(yīng)塑性自由度被激活)。

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

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

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

2.3 單元控制方程

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

圖9 建立的兩基本結(jié)點(diǎn)單元可作為模型的基礎(chǔ)單元，包含兩個(gè)剪切子單元及一個(gè)界面子單元，建立其虛功方程如下：

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

此外，考慮單元中各彈簧的非線性本構(gòu)關(guān)系可建立補(bǔ)充方程，對(duì)于進(jìn)入非線性狀態(tài)的單元內(nèi)力使用下式求解：

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

式中：ke為單元的初始彈性剛度矩陣；k′為描述單元結(jié)點(diǎn)力與非線性變形之間聯(lián)系的系數(shù)矩陣；為描述單元非線性行為的系數(shù)矩陣。具體通過下式計(jì)算：

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

3 結(jié)構(gòu)控制方程求解

將基礎(chǔ)單元控制方程進(jìn)行集成，可得整體結(jié)構(gòu)的隔離非線性控制方程：

式中： ΔX和 ΔF分別為整體結(jié)構(gòu)的結(jié)點(diǎn)位移增量向量和結(jié)點(diǎn)荷載增量向量；為整體結(jié)構(gòu)的非線性變形增量向量，其中僅包含了結(jié)構(gòu)中呈激活狀態(tài)的m個(gè)塑性自由度；Ke為整體結(jié)構(gòu)的初始彈性剛度矩陣，其規(guī)模為n×n階，階數(shù)等于結(jié)構(gòu)位移自由度數(shù)n；K′為描述整體結(jié)構(gòu)結(jié)點(diǎn)力與非線性變形之間聯(lián)系的系數(shù)矩陣，其規(guī)模為n×m階；為描述整體結(jié)構(gòu)材料非線性行為的系數(shù)矩陣，其規(guī)模為m×m階。當(dāng)結(jié)構(gòu)發(fā)生局部材料非線性變形時(shí)，激活的塑性自由度數(shù)目m將遠(yuǎn)小于結(jié)構(gòu)自由度數(shù)目n(即m＜＜n)。

式(38)表明：當(dāng)m遠(yuǎn)小于n時(shí)，任意時(shí)刻結(jié)構(gòu)的切線剛度矩陣可表示為初始彈性剛度矩陣的低秩攝動(dòng)形式。

采用完全Newton-Raphson 迭代法進(jìn)行結(jié)構(gòu)非線性分析，引入Woodbury 公式(式(39))在每個(gè)迭代步中對(duì)相應(yīng)式(38)進(jìn)行高效求解：

式中，KSchur為Schur 補(bǔ)矩陣：

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

4 數(shù)值算例

4.1 試驗(yàn)驗(yàn)證

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

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

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

圖10 砌體墻片試驗(yàn)滯回曲線Fig. 10 Experimental behaviour of simple piers

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

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

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

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

4.2 動(dòng)力非線性分析

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

圖14 結(jié)構(gòu)平面簡(jiǎn)圖Fig. 14 Plane of the masonry structure

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

圖16 約束砌體結(jié)構(gòu)數(shù)值分析模型Fig. 16 Discrete macro-element model of the masonry structure

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

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

圖18 結(jié)構(gòu)基底剪力時(shí)程曲線Fig. 18 Base shear force curve of the masonry structure

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

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

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

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

圖21 時(shí)間復(fù)雜度時(shí)程曲線Fig. 21 Time-history curves of time complexity

5 結(jié)論

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

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

(2) 通過分離結(jié)構(gòu)彈塑性狀態(tài)使砌體結(jié)構(gòu)非線性分析時(shí)僅對(duì)小規(guī)模非線性矩陣進(jìn)行更新和分解，可以避免傳統(tǒng)變剛度法中對(duì)整體切向剛度矩陣的更新與分解，該方法的時(shí)間復(fù)雜度明顯低于傳統(tǒng)變剛度法，可以實(shí)現(xiàn)砌體結(jié)構(gòu)非線性問題的高效求解。