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

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

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

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

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

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

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

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

1.1 基本假設(shè)

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

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

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

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

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

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

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

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

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

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

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

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

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

記作：

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

1.2 相對位移分解

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

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

圖3 單元法向應(yīng)力-相對位移關(guān)系Fig.3 Normal stress-relative displacement of Goodman element

圖4 單元剪切應(yīng)力-相對位移關(guān)系Fig.4 Shear stress-relative displacement of Goodman element

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

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

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

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

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

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

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

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

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

1.3 單元控制方程

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

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

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

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

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

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

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

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

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

其中：

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

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

式中：

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

2 考慮土-結(jié)相互作用結(jié)構(gòu)地震反應(yīng)分析

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

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

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

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

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

其中：

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

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

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

3 數(shù)值驗證

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

表1 接觸面本構(gòu)參數(shù)Table 1 Constitutive parameters of interface

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

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

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

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

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

4 數(shù)值算例

4.1 算例概況

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

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

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

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

表3 混凝土材料參數(shù)Table 3 Material parameters of concrete

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

4.2 結(jié)果分析

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

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

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

圖15 法向相對位移反應(yīng)Fig.15 Normal relative displacement response

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

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

圖16 剪切相對位移反應(yīng)Fig.16 Shear relative displacement response

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

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

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

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

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

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

5 結(jié)論

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

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

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