基于各向異性損傷本構模型的砌體結構抗震分析

富秋實， 錢 江

（同濟大學 土木工程與防災國家重點實驗室，上海200092）

1 引 言

砌體結構是土木工程中最古老的結構形式之一，為了有效地保護歷史建筑遺產和人員生命安全，針對既有砌體結構的結構分析手段顯得十分重要。一般，砌體結構分析主要是使用基于經驗的分析方法。過去的幾十年研究者對砌體結構的數值模擬方法進行了大量研究。另一方面，由于砌體材料特殊的各向異性性質，采用其他領域（如混凝土、巖石和復合材料）的數值模擬方法存在一定困難。所以，針對砌體的各向異性本構模型的研究開發具有重要意義。

在各向異性材料本構模型的研究領域，除了直接提出各種各向異性材料的屈服準則函數以及損傷準則函數［1］外，還有一類策略是采用基于應力張量分量線性變換（linear transformation）的方法。使用線性變換的方法便于保持屈服準則函數和損傷準則函數等的外凸性，這對于在實踐中保持本構算法的數值穩定性十分重要。因此，這類方法得到了較為廣泛的關注。為了有效地利用已有損傷理論，開發能夠考慮砌體材料各向異性特性的損傷模型，本文采用基于應力張量線性變換的方法。

應力張量的線性變換方法最早由Sobotka［2］和Boehler等［3］提出。Barlat等［4，5］將各向同性屈服函數與線性變換后應力張量的主應力相結合，用于解決正交各向異性材料的平面應力問題，并將這一方法應用于三維應力狀態。Karafillis等［6］將該方法推廣為各向同性塑性等效IPE（isotropic plasticity equivalent）理論，將線性變換方法與更加一般化的屈服函數相結合，可用于其他類型的材料。Barlat等［7］使用該方法提出了基于應力偏量線性變換的各向異性屈服函數，用于模擬金屬板在軋制方向和橫向的各向異性塑性特性。

砌體材料各向異性本構模型的已有研究，除了直接提出各項異性屈服函數，也有應用張量線性變換方法的研究［8，9］。這些研究使用的方法稱為映射張量法（mapped tensor method）。

映射張量法最早由Betten［10］提出，用于解決蠕變相關問題。其假設存在一個真實的應力空間和一個虛擬的應力空間，兩個應力空間之間存在映射關系，由一個四階映射張量定義。存在于真實應力空間的各向異性材料的受力狀態可以映射到存在于虛擬應力空間的虛擬各向同性受力狀態。可以看到，Betten的映射張量法的思想本質上與前文介紹的線性變換方法基本一致。Oller等［11，12］對這一映射張量法進行了進一步的發展。

本文進一步完善了應力張量的線性變換方法，提出了適應塑性以及損傷演化程度的動態線性變換張量，并將其應用于砌體材料本構模型的開發。利用該本構模型對砌體結構的振動臺試驗進行了模擬。模擬結果驗證了本文提出的彈塑性損傷本構模型的有效性，以及在進行結構非線性分析方面的優越性。

2 動態變換張量

設四階動態線性變換張量為Aσ，（Vk）， ∈｛p，d｝。內變量Vk是描述材料內部不可逆變化進程的變量，如塑性強化參數等。上標σ和p的含義為用于塑性本構的應力轉換張量，σ和d的含義為用于損傷本構的應力轉換張量。不同于已有的Oller和 Pelà模型［9，11，12］，本文并沒有應變轉換張量。但為了與文獻的應變轉換張量相區別，下文將依然保留上標中的σ。

根據應力轉換張量的定義，有

式中σ＊為變換后的應力張量。設未知其具體表達的各向異性屈服／損傷面函數F為

設已知形式的各向同性本構的屈服函數和塑性勢函數表達式為

此時根據定義顯然有

即得到了未知函數F的隱式表達。另一方面，由于不同方向各向異性強化的特性記入變換張量Aσ，（Vk）中，因此隨著強化參數Vk的演化，在F＊呈現各向同性強化的同時，F可以呈現各向異性強化。

3 各向異性彈塑性損傷模型建立

3.1 塑性演化

在以往研究中，對于混凝土等摩擦性材料的屈服函數一般采用 Mohr-Coulomb或Drucker-Prager屈服函數，但與試驗結果并不是很接近。本節同時采用Rankine和Drucker-Prager形式的雙屈服準則，模擬砌體材料受拉及受壓屈服特性，基于動態線性變換張量，建立針對砌體材料的正交各向異性塑性修正方法。

（1）受拉Rankine形式屈服面

各向同性受拉屈服函數采用Rankine形式為

引入四階受拉塑性轉換張量Aσ，pt，則可由式（1～3）得到受拉Rankine形式屈服面的各向異性隱式表達。

（2）受壓Drucker-Prager形式屈服面

各向同性受壓屈服函數采用Drucker-Prager形式為

式中α為控制雙軸受壓強度特性的參數，單軸受壓狀態下有效應力閾值fc＊簡化為硬化參數的線性函數。

引入四階受壓塑性轉換張量Aσ，pc，可由式（1～3）得到受壓Drucker-Prager形式屈服面的各向異性隱式表達。

3.2 損傷演化

為了反映在受拉和受壓狀態下材料的力學性能差異，損傷本構模型需要將應力張量σ［13］或者有效應力張量σ［14］分解為正負張量兩部分。在基于應變的有限元分析中，應力張量σ在迭代收斂前是未知量，基于應力張量的分解會極大地增加計算量和數值計算難度。另一方面，由于有效應力張量σ和彈性應變張量εe線性相關，基于有效應力張量的分解更利于建立應變的顯式表達。因此，本文采用有效應力張量的分解方法，建立損傷演化規則：

此外，損傷修正部分只需考慮彈性損傷力學模型，更新應力張量和損傷變量，而不需要考慮塑性變形的影響。因此，為了方便書寫，下文的理論推導部分將省略彈性應變張量εe的上標e。

（1）受拉損傷演化

受拉狀態下的各向同性損傷面函數采用Rankine形式：

τ＋＊具有應變張量的量綱，稱為受拉等效應力。為方便計算，損傷各向同性受拉準則函數寫為

r＋＊為受拉各向同性損傷面函數對應的等效應力閾值。

（2）受壓損傷演化

受壓狀態下的各向同性損傷面函數采用修正的八面體應力組合：

式中 參數K控制雙軸受壓狀態有效應力損傷閾值與各方向單軸受壓狀態下有效應力損傷閾值的比值，其取值由試驗確定。

（3）II型破壞損傷演化

砌塊-砂漿界面的II型破壞尚難以用損傷變量d＋和d－來描述，其滑動破壞機制也與受拉及受壓開裂有著明顯的區別。根據內變量理論的概念，不同機制的不可逆過程應該由不同的內變量來表示。因此，本節在前文理論的基礎上，直接引入一個損傷變量dII，來描述砌體材料灰縫的II型破壞。

首先，引入II型破壞等效應力：

式中σ12為沿灰縫滑動的剪應力，σ22為垂直于水平灰縫的正應力，II為II型破壞對應的摩擦角。相應的損傷準則函數為

式中r0，II為II型破壞等效應力的初始閾值，由剪切試驗確定。假定II型破壞損傷變量的顯式表達式為指數函數形式：

4 對比驗證

同濟大學于2011年對某砌體結構進行了振動臺試驗，該試驗原型結構為上海市某臨近地鐵的住宅［15，16］。由于振動臺試驗的能力限制以及砌體結構模型試驗的特殊性，試驗采用局部代表性結構來模擬原型結構，從實際結構中取出三層兩個房間作為代表性結構進行研究。由于代表性結構的質量（690t）依然遠超出振動臺的能力范圍，所以本試驗采用代表性結構1／4的縮尺模型進行振動臺試驗研究。

地震波加載方向為模型的橫墻方向，即X向。模型砌體部分材料采用小磚塊和砂漿，構造柱混凝土采用砂漿，鋼筋采用鍍鋅鐵絲，原型結構樓屋面厚度為110mm，1／4縮尺后樓屋面板厚為27.5mm，考慮到施加質量塊時固定砂漿會對樓板有一定程度的加強，并且樓板并不是試驗模擬的重點，為使樓板能夠承受附加質量引起的合重，樓屋面板厚取40mm，混凝土部分用砂漿，鋼筋采用鍍鋅鐵絲進行模擬。砂漿澆注時用小振搗器振搗和小鐵釬插搗，保證模型澆筑質量。模擬砌體部分的小磚塊由MU10磚塊在三個方向上一分為四，即每塊磚切成64塊模型小磚塊，砂漿用M2.5混合砂漿。

4.1 計算模型概況

通過商業有限元程序ANSYS的Usermat子程序，將前文所建立的各向異性本構模型嵌入ANSYS程序中，用來進行砌體結構非線性分析。程序編寫基于平面應力狀態，可用于平面應力單元以及殼單元（平面應力＋板彎曲）。

利用ANSYS程序進行三維空間建模。砌體墻、樓板及模型底座采用SHELL單元模擬，混凝土梁和柱采用BEAM單元模擬。計算模型共計單元數3392，節點數3494。計算采用的主要材料參數為x向彈性模量E1＝1806MPa，抗拉強度＝0.6MPa，抗壓強度＝4.8MPa；y向彈性模量E2＝1508MPa，抗拉強度＝0.4MPa，抗壓強度f2－＝4.0MPa。計算模型如圖1所示。

4.2 8度罕遇烈度El-Centro波作用下結構響應結果對比

為了盡可能使計算工況與試驗工況保持一致，計算時的加速度輸入采用相應工況的試驗臺面記錄。

（1）樓層加速度對比

表1及圖2為8度罕遇El-Centro波下，加速度試驗與數值分析包絡結果對比。

圖3～圖6為8度罕遇El-Centro波下，加速度試驗與數值分析加速度時程結果對比。

圖1 有限元模型Fig.1 FEM model

表1 加速度試驗與數值分析包絡結果對比Tab.1 Comparison of acceleration results obtained from test and simulation

可以看出，除了模型底板處的動力響應結果與試驗偏差略大以外，整體結構的相應幅值與試驗結果比較接近；相位結果有一定偏差，其原因可能為計算模型未考慮混凝土底梁的非線性行為，但在試驗過程中混凝土底梁進入了非線性狀態。

圖2 加速度試驗與數值分析包絡結果對比Fig.2 Comparison of acceleration results obtained from test and simulation

圖3 一層加速度反應時程結果對比Fig.3 Comparison of acceleration time history curves of the first floor

圖4 二層加速度反應時程結果對比Fig.4 Comparison of acceleration time history curves of the second floor

圖5 三層加速度反應時程結果對比Fig.5 Comparison of acceleration time history curves of the third floor

圖6 頂板層加速度反應時程結果對比Fig.6 Comparison of acceleration time history curves of the roof

（2）樓層相對位移對比

表2和圖7為8度罕遇El-Centro波作用下，樓層相對位移試驗與數值分析包絡結果對比。兩者結果除了2層處相差較遠外，其他位置結果比較接近。考慮本模型在2層處沒有附加質量和剛度突變，計算所得的位移包絡曲線形狀更加符合結構概念，因此，認為計算得到的結果具有較好的可信度。

（3）損傷分布情況對比

計算結構損傷分布情況如圖8所示，損傷主要分布在模型的底層。試驗結果顯示，在8度罕遇El-Centro波作用下，試驗模型主要破壞模式為橫墻1／2高度處起始的斜向階梯型裂縫（試驗過程中能聽到明顯的開裂聲音），相應的試驗模型裂縫照片和計算模型損傷分布對比情況如圖9所示，可見計算模型的損傷分布與試驗模型裂縫開展情況基本相符。同時，縱墻的墻體底部也發生了較為輕微的水平開裂，試驗模型裂縫照片和計算模型損傷分布對比情況如圖10所示，兩者結果也基本相符。

表2 樓層相對位移試驗與數值分析包絡結果對比Tab.2 Comparison of maximum story relative displacements

圖7 樓層相對位移試驗與數值分析包絡結果對比Fig.7 Comparison of maximum story relative displacements

圖8 計算模型整體損傷分布情況Fig.8 Damage distribution of FEM model

圖9 底層橫墻裂縫／損傷分布對比Fig.9 Crack／damage distribution of bottom lateral wall

5 結 論

本文進一步完善了應力張量的線性變換方法，提出了適應塑性以及損傷演化程度的動態線性變換張量，建立了以該變換張量為基礎發展各向異性彈塑性本構模型以及損傷本構模型的一般方法，并將其應用于砌體材料本構模型的開發。

基于開發的砌體各向異性彈塑性損傷本構模型，本文對砌體結構的振動臺試驗進行了模擬。模擬結果進一步驗證了本文提出的彈塑性損傷本構模型的有效性，以及在進行結構非線性分析方面的優越性，即除了能夠獲取結構的反應外，還能實時提供結構中的損傷分布狀態信息，找出結構的薄弱部位，并據此合理地設計結構或進行相應的結構修復。