考慮土-結相互作用的大型結構高效地震分析方法

胡靜靜，余丁浩，李 鋼，王 睿，張 晗，蘇 璞

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

隨著經濟社會的迅速發展和科學技術的不斷突破，近年來我國逐漸涌現大量的大型復雜工程結構，如大型城市綜合體、地下核電站、超高層建筑物等。這些結構具有如下特點[1]：① 體量大且構成復雜；② 土與結構接觸部位剛度突變明顯、相互作用效應顯著且不容忽略。這些大型工程結構在地震下的安全問題關系著國計民生[2]，因此對此類結構進行地震反應分析具有重大意義，特別是開展考慮土-大型結構相互作用的高效非線性分析方法研究顯得尤為重要。

目前，考慮土-結動力相互作用的數值分析方法主要有直接法和子結構法[3]。子結構法是先把地基和上部結構分別離散分析，再根據接觸界面上的連續性條件考慮兩個子結構間的相互作用，該方法的優點是靈活性大、計算量小，國內外相關學者已基于此方法進行了大量研究[4-5]。然而，該方法前期主要限于線性或者等效線性問題求解，為將其擴展到非線性地震反應分析方面，WOLF[6]總結了在時域中進行土-結相互作用非線性分析的步驟，提出線性和非線性體系均適用的頻域-時域混合分析程序，但對于非線性體系來說該方法計算結果僅是一種有限近似；傅敏紅等[7]基于子結構法提出了數值解和解析解耦合的土-結相互作用求解方法，可將其拓展到非線性計算領域，但該方法求解過程中基礎和地基土必須分別假定為剛性和彈性，因而其適用性有限。由于子結構法是在疊加原理的基礎上建立的，僅可用于某些特定條件下的非線性地震反應分析，這使得此類方法的應用受到較大限制[8]。直接法也稱整體有限元法，該方法將結構與較大范圍的近場地基土域看作整體進行建模分析，可以更完整、真實的考慮結構在地震作用下的動力響應，能夠很好的模擬土與結構界面的分離、滑移等強動力非線性接觸行為。整體有限元法[9]已成為處理土-結相互作用問題的最有效方法，該方法最關鍵的環節是如何處理土與地下結構之間接觸面的動力學行為。地震作用下，土與周圍結構的變形不一致，土-結接觸面上發生的相對位移、張開、重新閉合等強接觸非線性行為[10]直接影響接觸面附近土體與結構的動力響應。目前有限元法在處理土與結構接觸面的動力接觸行為時主要有兩種方法[11]：一是接觸力學法，即通過直接定義不同介質之間接觸表面對的力學傳遞特性，建立接觸面力傳遞的力學模型和接觸方程，通過接觸算法求解接觸方程，代表性求解算法有Lagrange 乘子法[12-13]、Penalty法[14-15]等；另一種是接觸面單元法，即通過在兩種接觸的介質之間建立接觸面單元，同時賦予接觸面單元特殊的應力-應變關系來模擬接觸面的力學行為，代表性方法有Goodman 單元[16]和薄層單元[17]等。隨著有限元技術的不斷突破，接觸面單元法憑借其概念簡單、能夠與普通單元一樣進行程序化實現等優勢已逐漸發展為土-結接觸面數值分析的主流方法[18]。

隨著工程結構規模的擴大，采用整體有限元法求解考慮土-結相互作用的大型結構非線性地震反應時，由于需同時對主體結構、大范圍土域和接觸面進行建模，其單元數與結點自由度數通常非常巨大，這就對計算效率提出了更高要求。目前，結構非線性分析主要采用變剛度求解方法，該方法需要在迭代求解過程中對結構整體剛度矩陣進行實時更新和分解，對于大型結構而言，大規模剛度矩陣的時變特性使得非線性分析計算效率較低。為提高非線性分析效率，眾多學者分別從不同角度提出了多種加速計算方法，如方程組高效求解方法、快速迭代方法、并行計算方法等，這些方法或者通過提高大規模線性方程組求解速度來改善計算效率，或者通過減少剛度矩陣更新次數來提高效率，或者通過降低方程組規模來提高計算效率。地震作用下結構通常并非全部進入非線性狀態，而是僅在某些局部區域出現非線性，而結構的其他大部分區域則始終保持為彈性狀態，眾多學者基于工程結構的局部非線性特征提出了多種高效分析方法，主要有多尺度方法[19-21]、子結構方法[22-24]、重分析方法[25-26]、擬力法[27-29]等，然而，這些方法適用性通常受限，難以直接用來進行考慮土-結構相互作用的大型結構地震反應分析。近年來，通過將擬力法中的變形分解思想進行拓展，李鋼等提出了一種求解局部材料非線性問題的新型高效通用有限元算法，即隔離非線性有限元法[30-33]，該方法通過對材料應變進行分解并在單元內構造非線性應變場函數，能夠在非線性迭代求解過程中保持結構彈性剛度矩陣不變，避免了傳統方法由于結構整體剛度矩陣時變[34-35]導致的計算效率低下，在保證精細化模擬的同時大幅提高計算效率，目前已開發出基于該方法的多種單元模型，如纖維梁單元[36]、殼單元[37]、實體單元[38]等，但還未涉及土-結構相互作用問題的高效非線性分析。雖然土-結構接觸區域在地震作用下可能出現較大范圍的強接觸非線性變形，但此區域通常僅占整體計算模型的小部分，因此，考慮土-結構相互作用的非線性分析問題具有顯著的局部非線性特征，結合適當的局部非線性算法有助于顯著提升此類問題的分析效率。

本文提出一種考慮土-結相互作用的大型結構高效地震反應分析方法，該方法以隔離非線性有限元法為理論基礎，采用無厚度三維Goodman 單元模擬土-結接觸行為，通過將單元相對位移分解為初始彈性和非線性兩部分，并引入額外非線性自由度建立單元非線性相對位移場，建立了能夠模擬土-結構相互作用的隔離非線性接觸面單元控制方程，隨后對主體結構、土體和土-結接觸模型的控制方程進行集成建立了隔離非線性的整體結構動力控制方程，最后引入Woodbury 公式求解，在非線性分析時僅對代表結構局部非線性行為的小規模Schur 補矩陣進行不斷更新和分解，大大提高了結構非線性分析效率，數值算例結果驗證了本文方法的有效性和高效性。

1 隔離非線性的土-結相互作用計算模型

為實現土-結構接觸行為的高效模擬，本節基于變形分解思想建立隔離非線性的接觸面單元模型，考慮到Goodman 單元簡單易于實現，并且與界面本構模型結合可較好地模擬土-結構接觸區域的滑移、張開和壓密等接觸非線性行為，目前在工程中得到廣泛應用[39-43]，本文以該單元為基礎對土-結構的接觸非線性行為進行描述。

1.1 基本假設

八結點Goodman 接觸面單元是由兩個四結點平面組成，如圖1 所示面Ⅰ和面Ⅱ，單元厚度e=0，圖中ξηζ 代表單元等參坐標系，兩接觸面之間任意一接觸點對(點1 和點2，以下簡稱點對)的受力變形關系可認為通過一組剪切方向和法向彈簧進行模擬，通過對這些彈簧定義適當的應力-相對位移本構關系即可實現接觸面間水平相對錯動、法向張開和閉合等行為的計算分析。

圖1 Goodman 單元工作原理Fig.1 Principle of the Goodman element

如圖2(a)所示，空間八結點無厚度Goodman接觸面單元的兩個接觸面對應結點的坐標完全相同，圖中xyz代表整體坐標系，即：

圖2 三維Goodman 單元Fig.2 Three-dimensional Goodman element

式中：xi、yi、zi是接觸面Ⅰ (結構表面)的單元結點i(i=1, 2, 3, 4)在整體坐標系內的坐標值；xi+4、yi+4、zi+4是接觸面Ⅱ (土體表面)的單元結點i+4(i=1, 2, 3, 4)在整體坐標系內的坐標值。平面問題的Goodman 單元本質上是一維單元[44]，因而本文中三維Goodman 單元本質上可以看作某種特殊的二維單元，基于此可將空間八結點接觸面單元進行等參處理，如圖2(b)所示，其中ξηζ 代表等參坐標系，等參坐標系下單元形函數可表示為：

式中，ξi和ηi分別為結點i(i=1, 2, 3, 4)對應的等參坐標值，即ξi、ηi=±1。

為確定整體坐標系xyz與等參坐標系ξηζ 之間的對應關系，將(ξ,η,ζ)空間內形狀規則的單元映射為(x,y,z)空間的形狀扭曲的單元，需引入如下等參變換方程：

八結點接觸面單元每個結點具有3 個自由度，在迭代過程中的任一線性化增量求解步中，單元的結點位移增量(整體坐標系下)可表示為：

單元兩片接觸面之間任意點對的相對位移增量為Δu=[ΔuΔvΔw]T，其可由單元變形矩陣B和結點位移增量Δq表示：

式中，I為3×3 的單位矩陣。

在接觸面單元中，需考慮沿其剪切方向和法向的相對位移和應力，取接觸面單元上任一點沿接觸面的法線方向為ye軸，如圖2(b)所示，建立局部坐標系xeyeze，xeze即為過該點的切平面[45]。在局部坐標系下，單元相對位移增量表示為Δue=[ΔueΔveΔwe]T，其中Δue和Δve對應局部xe向和ze向的剪切相對位移分量，用來描述土-結接觸面之間的相對錯動變形；Δwe對應局部ye向的法向相對位移分量，用來描述土-結接觸面之間的張開、閉合等變形。局部坐標系下的Δue可由如下表達式計算：

式中，L為坐標轉換矩陣。令LB=B'，則有：

在局部坐標系xeyeze下，設ksη和ksξ為當前計算步下單元中任意點對在xe向和ze向(剪切方向)的切線剛度系數，kζ為該點對ye向(法向)的切線剛度系數，則局部坐標系下該點的增量應力-相對位移關系可表示為：

記作：

式中：Δσi=[ΔτηΔτξΔσζ]T為應力增量向量，其中，Δτη和Δτξ分別為局部xe向和ze向的剪應力分量，Δσζ為局部ye向的正應力分量；Dt為該點對在當前計算步下的瞬時切線剛度矩陣。

1.2 相對位移分解

在任一線性化增量迭代計算步中，接觸面單元中任意點對的相對位移增量向量Δue可根據該點的初始彈性剛度分解為初始彈性和非線性兩部分：

相對位移增量向量中包含2 個剪切相對位移分量Δue、Δve和法向相對位移分量Δwe(即Δue=[ΔueΔveΔwe])。單元兩個剪切方向的行為并無耦合，但法向行為會對剪切方向行為產生影響。圖3 和圖4分別為單元法向和剪切方向的典型應力-相對位移關系圖，當兩接觸面閉合受壓時(法向相對位移we＜0)，單元法向行為可通過一較大的剛度kp進行模擬(如圖3 所示)，此時單元存在剪切摩擦行為，相應的剪切應力-相對位移關系如圖4 所示；當兩接觸面打開時(法向相對位移we＞0)，單元法向正應力和剪切應力均為0 值。

圖3 單元法向應力-相對位移關系Fig.3 Normal stress-relative displacement of Goodman element

圖4 單元剪切應力-相對位移關系Fig.4 Shear stress-relative displacement of Goodman element

圖5 剪切相對位移分解Fig.5 Decomposition of shear relative displacement

由于法線方向上的應力-相對位移關系在原點附近存在拐點，導致單元內任意點對在不同加載方向上初始剛度不唯一，不同初始接觸狀態下單元變形分解的具體實現方式存在差異，相應的，對于接觸面非線性狀態的判定標準也不相同。為此，本文建立了不同初始接觸狀態下的單元位移分解方法，具體論述如下：

1) 當單元中某點對的法向初始狀態為受壓，且在當前計算步下其法向依然受壓，但兩個剪切方向的瞬時切線剛度偏離其初始剛度，則認為該點對在剪切方向進入非線性狀態，即在剪切方向產生非零的非線性相對位移，但法向依然處于初始彈性狀態，即法向的非線性相對位移分量為0，此時單元非線性相對位移增量表示為：

應當指出的是，本文方法中初始剛度主要用來作為變形分解的依據，從而使得在計算求解過程中整體剛度保持不變。當考慮卸載剛度退化或強度下降時，本文方法依然適用，此時，每個增量步中的變形增量均可以初始彈性剛度為基準分解為兩部分，分解出的非線性變形將為非零值，對應的非線性自由度將依然處于激活狀態，并參與控制方程和相應求解公式的構造與計算。

2) 當單元中某個點對的法向初始狀態為受壓，且在當前計算步下其法向變為受拉，則此時該點對三個方向的切線剛度均偏離其初始剛度(均為0)，因而均會判定為進入非線性狀態，相應非線性相對位移增量的三個分量均為非零值。圖6給出了此時法向相對位移增量的分解示意圖，假設當前增量計算步下其應力狀態從a點變化到b點，則對應的非線性相對位移增量計算表達式為：

圖6 法向相對位移分解Fig.6 Decomposition of normal relative displacement

3) 當單元中某點對的法向初始狀態為受拉，且在當前計算步下其法向依然受拉，則此時該點對的三個方向均處為零剛度狀態，與初始狀態相同，因而均認為處于初始彈性狀態，相應非線性相對位移增量向量中的三個分量均為0，即=[0 0 0]。

4) 當單元中某點對的法向初始狀態為受拉，且在當前計算步下其法向變為受壓，則該點對三個方向的瞬時切線剛度均不再是初始狀態時的零剛度，因而均認為進入非線性狀態，相應各分量的非線性相對位移計算表達式的形式與式(13)和式(14)相同，為防止初始零剛度引起的計算奇異問題，本文在計算時使用一個接近于0 的小數近似代替各方向初始時刻的剛度值。

本文中對于非線性相對位移的定義是相對于初始狀態而言的，代表了單元任意點對當前狀態與其初始狀態的數值偏離程度，其作用主要是為了能夠始終使用恒定的初始剛度對單元任意時刻的應力狀態進行描述，以便于下文推導出能夠使結構彈性剛度矩陣在分析過程中始終保持不變的隔離非線性控制方程。

1.3 單元控制方程

任意狀態下單元的應力增量場σin可以用單元彈性相對位移增量場與彈性矩陣De的乘積來計算，結合式(10)可將其進一步表達為：

相較于傳統有限元方法，隔離非線性方法需在單元中額外設置若干個非線性變形插值點，以便于對單元非線性變形場進行模擬，本文以單元高斯點作為該單元的非線性相對位移插值點，并在每個接觸單元中設置4 個高斯積分點，如圖7所示，則單元的非線性相對位移場可通過如下插值形式給出：

圖7 單元非線性變形插值點Fig.7 Inelastic deformation interpolation point of the element

式中，I3為3×3 的單位矩陣。

將式(7)和式(16)代入式(15)可將任意狀態下單元的應力增量場Δσin進一步表示為：

在計算分析過程中，對于彈性單元，其相應的應力增量可表示為單元初始剛度矩陣與單元相對位移的乘積，而對于進入非線性狀態的單元，其任意點對在任意線性化增量迭代步中的應力增量Δσin又可通過材料瞬時切線本構矩陣Dt與總相對位移增量Δue乘積的形式求得(如式(9)所示)，結合式(9)和式(10)可得出如下表達式：

將式(11)與式(16)代入式(19)，整理后可將單元應力增量場的計算表達式重寫為如下形式：

引入虛功原理建立單元控制方程，其中本文接觸面單元的虛功方程為：

根據矩陣運算相關規則，將式(23)進行整理，即可建立接觸面單元的控制方程：

其中：

式中：ke為單元的初始彈性剛度矩陣，在計算過程中始終不變，其具體數值與單元初始狀態有關；k′和為與單元非線性有關的系數矩陣，可通過高斯積分方法進行數值求解。由于本文令非線性變形插值點與高斯積分點重合，可以證明插值函數矩陣C滿足Kronecker delta 性質[46]，因此將單元的第j個高斯積分點坐標代入非線性變形插值函數Ci后可得：

式中：Dt,i代表當前計算步下單元非線性插值點處的材料切線本構矩陣；Hi為第i個插值點處的積分權系數；為第i個插值點處的變形矩陣；|J|i為插值點處的雅可比矩陣行列式[45]，其計算表達式為：

式中：

在實際計算過程中，可能存在某些插值點的全部或個別非線性變形分量被判定處于初始彈性狀態的情況，如1.2 節相對位移分解情況1)和情況3)所述，此時這些插值點中相應非線性變形分量為0，在單元控制方程推導時應將這些變形分量直接剔除，在實際程序實現時，只需要將單元控制方程式(27)中相應非線性自由度及系數矩陣k′和中對應行、列直接刪除即可，這一過程可由程序根據單元狀態自動完成，當整個單元都處于線彈性狀態時，單元控制方程將變為keΔq=ΔFe。

2 考慮土-結相互作用結構地震反應分析

本文采用整體式模擬方法建立考慮土-結構相互作用的結構分析模型，以圖8 所示框架剪力墻高層結構為例進行說明，計算模型中包含主體結構、土體和土-結構接觸區域三部分，其中，主體結構的梁柱構件采用隔離非線性纖維梁單元[47]進行模擬；主體結構中的剪力墻采用隔離非線性殼單元[37]進行模擬；土體部分采用隔離非線性實體單元模型[48]進行模擬，該模型將材料應變分解為線彈性應變與非線性應變兩部分，并通過插值形式構造非線性應變場，可實現土體部分局部非線性行為的隔離表達和高效分析；土-結構接觸區域采用本文前述建立的隔離非線性接觸單元進行模擬。通過將不同單元的控制方程進行集成，可得如下形式的整體計算模型控制方程：

圖8 隔離非線性單元模型Fig.8 Element model of inelasticity-separated method

式中：Ke為結構整體初始彈性剛度矩陣，在計算過程中保持恒定不變；K'和由各單元控制方程中的系數矩陣k'和集成而來；ΔQ為結構的結點位移增量；ΔF為結構的荷載增量；集合了整個結構中所有非線性單元中的非線性插值點的非零非線性變形。K'與材料的非線性程度無關，僅與產生非零增量非線性變形的區域位置有關，矩陣不僅與產生非零非線性變形的區域位置有關，還與材料的非線性本構模型相關。結構控制方程表現為整體彈性剛度矩陣Ke和代表局部非線性變形的剛度項K'和相互分離的特征，設結構整體自由度數為n、非線性自由度數為m，則矩陣Ke、K'和的規模分別為n×n、n×m、m×m，由于控制方程中非線性相關系數矩陣在集成時僅考慮了局部區域內進入非線性的相關單元，m代表了結構非線性區域的規模，其值通常遠小于n(即m＜＜n)，這表明控制方程式(31)中非線性自由度和相關系數矩陣的規模通常極小。

對于考慮土體-結構相互作用的整體結構動力非線性分析，可列出增量形式的運動微分方程為：

式中：M為結構質量矩陣；C為結構阻尼矩陣；ΔQ¨(t) 、 ΔQ˙(t)分別為結構各時刻的相對加速度以及速度； ΔPg(t)為地震作用力；ΔF(t)為結構的恢復力向量。將式(31)代入式(32)，并結合Newmarkβ 法求解運動微分方程，可基于隔離非線性法建立考慮土體-結構相互作用的整體結構動力分析控制方程：

其中：

算法實現流程如圖9 所示?？梢钥闯?，在每次迭代時的單元狀態確定過程中，需對各單元的非線性狀態進行判定，這可通過將單元每個非線性插值點處的瞬時切線剛度與其初始彈性剛度對比實現，對于某個單元，若全部插值點處的瞬時切線剛度與初始彈性剛度相等，就表明該單元在當前增量計算步下處于初始彈性狀態，否則就代表該單元在當前計算步下處于非線性狀態。

圖9 算法流程圖Fig.9 Algorithm flow chart

此外，在每次迭代時，根據確定的非線性單元位置和數量，可直接確定非線性自由度數目，從而確定求解方程中非線性相關系數矩陣的規模及其中的非零元位置和數量，本文方法采用Fortran 編程，利用動態內存分配機制實時開辟和銷毀相關系數矩陣所需存儲空間?？梢钥闯?，相較于傳統方法，本文方法雖然計算效率更高，但需增加額外的存儲開銷，然而，由于計算求解過程中涉及到的非線性相關系數矩陣均具有高度的稀疏性[34]，且本文方法采用僅存儲非零元的行壓縮格式(compressed sparse row，CSR)進行矩陣存儲，因此額外增加的存儲開銷并不顯著。

3 數值驗證

為驗證本文方法準確性和高效性，建立如圖10所示的接觸摩擦算例模型，上部塊體作用有給定豎向荷載P1和單調遞增的水平荷載P2，考慮兩個計算工況，其中工況1 中P1荷載為50 kN，工況2中P1荷載為100 kN。上部塊體與下部塊體均采用實體單元模擬，彈性模量分別為2.1×1010Pa 和5×107Pa，泊松比均為0.3。上、下塊體接觸面采用本文建立的接觸面單元模擬，接觸本構采用Clough-Duncan 雙曲線模型[50]，具體本構模型參數取值詳見表1。

表1 接觸面本構參數Table 1 Constitutive parameters of interface

圖10 接觸摩擦塊模型 /mFig.10 Contact friction block model

采用本文方法對該模型進行分析，兩個工況下計算所得接觸面平均剪切應力-相對位移曲線如圖11 所示，圖中與理論解進行了對比，其中理論解為基于Clough-Duncan 模型計算出的接觸面極限剪切應力(即剪切相對位移u→∞時的剪切應力值)，表2 給出數值解(即本文方法加載至最大位移處對應的平均剪應力)與理論解的相對誤差。由圖11 和表2 可見，本文方法計算結果與理論解吻合良好。

表2 本文方法數值解與理論解對比Table 2 Comparison between numerical solution and theoretical solution

圖11 剪切應力-相對位移曲線Fig.11 Shear stress-relative displacement curve

為了驗證本文方法的高效性，利用有限元計算軟件ABAQUS 對該算例進行分析，模型規模相同條件下兩個工況計算時間分別為344 s 和342 s，而本文方法的計算時間分別為35 s 和47 s，僅為ABAQUS 的10.2%和13.7%。

4 數值算例

4.1 算例概況

采用本文方法對某大型城市綜合體結構進行考慮土體-結構相互作用的地震反應分析，結構類型為框架核心筒體系，結構示意圖如圖12 所示。該結構地上31 層、地下2 層，裙房地上部分共5 層，地上結構高度99.45 m，地下埋深-10.5 m，平面尺寸約為126 m×67.2 m，建筑總面積為115 506.72 m2。截取的地基土域平面尺寸為378 m×201.6 m，深度取60 m。對結構進行雙向地震動加載，為使分析過程中模型結構出現較為顯著的非線性變形，以便于對本文方法進行驗證，x方向和z方向的地震波峰值加速度分別為2.2 m/s2和1.87 m/s2，兩個方向的地震加速度時程如圖13 所示。

圖12 整體土-結有限元模型Fig.12 soil-structure finite element model

圖13 地震波加速度記錄Fig.13 Acceleration record of the earthquake

該綜合體的整體土-結計算分析模型如圖12 所示，其中主體結構的梁、柱構件采用隔離非線性纖維梁單元模擬，其單元數為80 658；剪力墻結構采用隔離非線性殼單元模擬，其單元數為76 391；土體采用隔離非線性實體單元模型模擬，其單元數為138 160；土-結相互作用采用本文隔離非線性接觸面單元模擬，其單元數為1116。整體模型的單元總數為296 325，總自由度數為615 513。模型中梁、柱構件的混凝土采用修改后的修正Kent-Park 本構模型[34]進行模擬，剪力墻混凝土材料采用修正斜壓場模型，混凝土材料參數見表3；鋼筋采用Dodd 和Cooke[51]提出的本構模型；地基土域的材料采用Drucker-Prager 模型[52-53]，鋼筋及土體相關材料參數見表4；土域邊界采用滾軸邊界進行處理[54]；土體與結構接觸面的本構關系選用Clough-Duncan 雙曲線模型[50]，其中界面摩擦角取值為31°，其他參數取值均同表1。本文算例模型對于樁基礎的模擬采用“復合式分析”方法，即將樁彌散于土體中，把群樁基礎視為樁、土合一的復合材料[55]，以降低建模和分析的復雜度。此外，本文算例中，主體結構的地下部分外側均為框架-剪力墻，其中梁柱構件采用隔離非線性纖維梁單元模擬，梁柱與剪力墻通過共結點方法整體建模，接觸單元僅連接墻單元和土體單元。

表3 混凝土材料參數Table 3 Material parameters of concrete

表4 鋼筋及土體材料參數Table 4 Material parameters of steel and soil

4.2 結果分析

對建立的考慮土-結相互作用的整體模型分別采用傳統非線性分析方法和本文方法進行地震反應分析，其中傳統分析方法通過在每次迭代時對整體剛度進行更新分解計算結構位移響應。

任選土體與結構接觸面上四個參考點對，分別為圖14 中所示A、B、C和D四個位置，這四個點對的法向相對位移時程響應如圖15 所示，圖中法向相對位移是指結構與土體之間的相對開度，接觸面張開為正，可以看出，四個點對的初始狀態均為受壓接觸。由圖15 可以看出，在前5 s內，這四個參考點處接觸點對的法向相對位移均在0 附近，說明結構與土體在該時間段內是接觸狀態；在5 s 之后，接觸點對的相對位移數值呈現“忽高忽低”的波動現象，且其最大值開始逐漸增大，說明這些接觸點對附近的土-結構接觸面表現出開、合交替出現的行為，并在10 s 左右時法向相對位移達到峰值，此外，由圖15(d)可以看出，點D處接觸點對的法向相對位移在9.5 s 之后始終大于0，說明該點附近的接觸面在9.5 s 之后始終處于脫開狀態。

圖14 接觸面參考點對Fig.14 Reference point pair of interface

圖15 法向相對位移反應Fig.15 Normal relative displacement response

圖16 給出這四個點對的剪切相對位移時程曲線。綜合各點對的法向和剪切相對位移時程曲線，可知接觸面發生了法向的脫開和重新閉合、剪切方向的相對摩擦和錯位滑移等非線性行為，這也驗證了本文方法模擬接觸面變形的有效性。此外，圖15 和圖16 表明本文方法與傳統方法的分析結果基本一致，表5 給出本文方法和傳統方法的峰值點相對位移數據及其相對誤差值，可以看出，本文方法具有較高的計算精度。

表5 相對位移相對誤差Table 5 Relative error of relative displacement

圖16 剪切相對位移反應Fig.16 Shear relative displacement response

圖17 和圖18 分別給出了結構在x和z方向的結構頂點位移時程曲線和層間位移角曲線，為分析土-結構相互作用對結構動力響應的影響，建立不考慮土-結相互作用的計算模型并進行分析(圖中模型a 為圖12(a)所示主體結構模型，模型b 為圖12(c)所示土-結相互作用整體模型)，從圖17 和圖18 可以看出，采用本文方法和傳統方法對模型b 分析的結果基本一致。此外，由圖17(a)和圖17(b)可以看出，在第10 s～12 s 內兩種模型分析得到的結構頂點位移達到峰值，且在x方向模型a 的峰值位移小于模型b；由圖18(a)可以看出，模型a 除第7 層～第13 層的x方向層間位移角大于模型b 外，其他樓層層間位移角均小于模型b；由圖18(b)可以看出，模型a 的層間位移角均小于模型b。上述結果表明：不考慮土-結相互作用會低估結構的位移反應。

圖17 頂點位移時程曲線Fig.17 Time-history curves of top displacement

圖18 層間位移角時程曲線Fig.18 Time-history curves of inter-story displacement angle

圖19 給出了本文方法在整體土-結相互作用模型計算過程中的非線性自由度變化過程。從圖中可以看出，結構在前200 分析步內已進入了非線性，但非線性自由度數目很小。整個分析過程中產生的非線性自由度最大值為135 206，占總自由度數(615 513)的22%，而每個增量步的平均非線性自由度為59 643，僅占總自由度數(615 513)的9.7%。這說明，本文方法在進行該結構地震反應分析過程中需要實時更新和分解的Schur 補矩陣階數遠小于結構整體剛度矩陣階數，Woodbury 公式的使用可以顯著提升結構分析效率。

圖19 非線性自由度時程曲線Fig.19 Time-history curves of inelastic degrees of freedom

本文方法的高效性可通過與傳統非線性分析方法的計算耗時進行對比來驗證。其中，傳統分析方法需在每個迭代步進行一次整體剛度矩陣的更新和分析，兩種方法計算過程均采用相同平臺，計算采用的硬件設備為小型服務器(CPU 為2 個Intel Xeon Gold 6254，內存為384 GB)。采用本文方法對該算例進行分析共耗時約6311 s；而采用傳統非線性分析方法則需耗費約42 072 s，為本文方法所需計算時間的6.7 倍。

5 結論

為實現大型結構考慮土體-結構相互作用的高效地震反應分析，本文建立了基于隔離非線性的土-結構接觸行為非線性計算模型：首先根據初始接觸狀態建立了接觸面單元的相對位移分解方法；隨后引入額外非線性自由度并結合Goodman 接觸單元格式對接觸面上的非線性變形進行模擬，能夠保證單元剛度矩陣在分析過程中恒定不變；進一步采用Woodbury 公式進行求解，有效提升了此類問題的分析效率。本文得到結論如下：

(1) 隔離非線性法要求根據初始彈性剛度進行變形分解，以便于保持結構剛度的恒常性及實現局部非線性行為的隔離表達，然而，單元法向應力-相對位移關系在零點附近存在拐點，導致單元內任意點對在不同加載方向上的初始剛度不唯一，為在分析過程中始終使用恒定的初始剛度對單元任意時刻的應力狀態進行描述，本文分別對不同初始狀態下的單元變形分解的非線性相對位移進行了定義。數值算例分析結果表明，本文建立的三維接觸面單元模型可對土與結構接觸面處的非線性行為進行有效模擬，且不考慮土-結相互作用會低估上部結構的位移響應。

(2) 計算結果顯示，分析過程中出現的非線性自由度較少，采用Woodbury 公式對結構控制方程進行求解時能夠避免傳統非線性分析方法在迭代求解過程中所需的大規模整體剛度矩陣的實時更新和分解，僅需對一個代表局部非線性行為的小規模Schur 補矩陣進行迭代更新，在保證精度的同時顯著降低求解時間，對大型結構分析問題具有較高的效率優勢。