第一作者姜忻良男，博士，教授，博士生導師，1951年生

郵箱：jiangxinliang@126.com

土-結(jié)構(gòu)非線性相互作用混合約束模態(tài)實施方法

姜忻良，張海順

(天津大學建筑工程學院濱海土木工程結(jié)構(gòu)與安全教育部重點實驗室，天津300072)

摘要：分析、探討線性-非線性混合約束模態(tài)綜合法用于土-結(jié)構(gòu)相互作用的動力分析實施方法。據(jù)結(jié)構(gòu)存在局部非線性特性，提出對整體結(jié)構(gòu)劃分為若干個線性、非線性子結(jié)構(gòu)，對線性子結(jié)構(gòu)只需一步提取特性矩陣并進行方程降階，而對非線性子結(jié)構(gòu)則逐步提取降階后的等效特性矩陣。通過組裝各子結(jié)構(gòu)達到對整體結(jié)構(gòu)特性矩陣的降階處理。對非線性子結(jié)構(gòu)憑借在小段時間Δt內(nèi)利用分段等效線性化手段，通過ANSYS二次開發(fā)工具UPFs編輯寫成可執(zhí)行程序文件，并與MATLAB聯(lián)立計算，實現(xiàn)特性矩陣提取。結(jié)果表明，該方法求得結(jié)構(gòu)的各階頻率和動態(tài)響應時程曲線與ANSYS直接計算吻合良好。可逐步提取等效彈塑性特性矩陣，為采用變化的特性矩陣對結(jié)構(gòu)進入塑性階段深層次分析研究提供有效方法。

關鍵詞：土-結(jié)構(gòu)相互作用；約束模態(tài)綜合法；勢能判據(jù)截斷準則；局部非線性；彈塑性剛度矩陣

收稿日期：2013-11-08修改稿收到日期：2014-03-31

中圖分類號:TU317.1文獻標志碼：A

基金項目：國家“863”計劃資助項目(2013AA040103)；河南省高等學校精密制造技術與工程重點學科開放實驗室開放基金資助(PMTE201308A)

Mixed constraint modal method for nonlinear soil-structure interaction

JIANGXin-liang，ZHANGHai-shun(School of Civil Engineering, Tianjin University/Key Laboratory of Coastal Civil Engineering Structure and Safety, Ministry of Education, Tianjin 300072, China)

Abstract:Here, the implementation technique of the linear-nonlinear mixed constraint modal synthesis method was analyzed for soil-structure dynamic interaction. According to the local nonlinear feature, the whole structure was divided into several linear and nonlinear sub-structures. Those linear sub-structures just required only one step to extract the characteristic matrices and reduce equation order, while for nonlinear sub-structures, they needed step by step to extract the reduced order equivalent characteristic matrices. By assembling the above sub-structures, the order of the whole structural characteristic matrix was reduced. For the nonlinear sub-structures, using the piecewise equivalent linearization method within a small time interval Δt, an executable program file to extract the characteristic matrices was edited by means of ANSYS secondary development tools UPFs and the simultaneous computation was done with MATLAB. Results show that the modal frequencies and dynamic response time history curves obtained with the above method agree well with those obtained with ANSYS direct calculation; using this method can extract the equivalent elastic-plastic feature matrices step by step, and it provides an effective tool for further studying the plastic phase dynamic analysis of structures with time-varying characteristic matrices.

Key words:soil-structure interaction; constrained modal synthesis method; potential energy mode cut-off criterion; local nonlinear; elastic-plastic stiffness matrix

解決地基土-結(jié)構(gòu)動力相互作用(SSI)問題時，有限元法為離散化強有力工具，但因其自由度太多致直接求解十分困難。尤其三維問題，結(jié)構(gòu)越復雜，自由度數(shù)越多，即便用人工邊界適當減小土體區(qū)域范圍，仍會耗費大量計算成本。因采用有限元技術進行動力分析時需進行結(jié)構(gòu)各單元特性矩陣總裝，導致整體矩陣龐大，一般微機難以滿足計算要求。因此可借助劃分動態(tài)子結(jié)構(gòu)方法，利用約束模態(tài)綜合法實現(xiàn)方程降階，為利用微機進行大型結(jié)構(gòu)有限元動力分析提供了可行、經(jīng)濟手段。而方程降階僅適用于線彈性結(jié)構(gòu)，對非線性彈塑性結(jié)構(gòu)則不適用。

本文基于動態(tài)子結(jié)構(gòu)法中約束模態(tài)綜合法(CMS)與勢能判據(jù)截斷準則，據(jù)土-結(jié)構(gòu)相互作用模型存在局部非線性特性，提出將整體結(jié)構(gòu)劃分為若干線性子結(jié)構(gòu)、非線性子結(jié)構(gòu)，對線性子結(jié)構(gòu)只需一步提取特性矩陣并進行方程降階，而對非線性子結(jié)構(gòu)則逐步提取其等效特征矩陣，通過組裝子結(jié)構(gòu)達到對整體結(jié)構(gòu)特性方程降階處理，繼而對地基土-結(jié)構(gòu)相互作用模型進行特征值方程計算及動力時程有限元分析，并將計算結(jié)果數(shù)值曲線與有限元直接法對比，分析所提方法的高效可行性。

1基本理論及方法闡述

動態(tài)子結(jié)構(gòu)法理論基礎為Rayleigh-Ritz法，是縮減自由度的有效方法，其將大型復雜結(jié)構(gòu)劃分為若干子結(jié)構(gòu)，分析研究每個子結(jié)構(gòu)的動力特性，并利用經(jīng)驗手段或理論方法確定其保留的低階主要模態(tài)信息，再據(jù)每個子結(jié)構(gòu)交界面的協(xié)調(diào)關系組裝成整體結(jié)構(gòu)的動力特性。此方法用于分析較少自由度的整體結(jié)構(gòu)，使大型復雜結(jié)構(gòu)整體動力特性問題得以解決。對采用何種經(jīng)驗手段或理論方法確定其保留的低階主要模態(tài)信息既能滿足工程需求又可靠性高問題，姜忻良等[1]在基于動態(tài)子結(jié)構(gòu)法中約束模態(tài)綜合法下提出勢能判據(jù)截斷準則進行子結(jié)構(gòu)主模態(tài)截斷，并推論驗證由各子結(jié)構(gòu)勢能范數(shù)隨截取主模態(tài)階數(shù)變化曲線判斷，當截取階數(shù)取合適值時，子結(jié)構(gòu)勢能將趨于收斂，即為最佳子結(jié)構(gòu)截取階數(shù)。

動態(tài)子結(jié)構(gòu)法與勢能判據(jù)截斷準的建立則均基于振型疊加法基礎，該類方法僅適合求解線性結(jié)構(gòu)系統(tǒng)動力問題，使動態(tài)子結(jié)構(gòu)法對非線性結(jié)構(gòu)體系進一步應用受到限制。研究表明，大量實際工程中整體結(jié)構(gòu)在外部荷載作用下并非全部構(gòu)件都進入非線性階段，而僅在某些區(qū)域才出現(xiàn)非線性特征，即結(jié)構(gòu)存在局部非線性特點。同樣土-結(jié)構(gòu)相互作用體系在地震動激勵下也存在局部非線性區(qū)域特點。在地震響應分析問題中，僅與上部結(jié)構(gòu)鄰近地基土區(qū)域會產(chǎn)生塑性應變，而遠離上部結(jié)構(gòu)的地基土區(qū)域在整個加載過程中始終處于線性階段。故考慮利用線性-非線性混合約束模態(tài)綜合法，將原不易進入非線性階段區(qū)域劃分成若干線性子結(jié)構(gòu)，再利用勢能判據(jù)截斷準則縮減自由度，而不必在整體結(jié)構(gòu)模型中反復迭代計算；而對非線性子結(jié)構(gòu)通過改變其各單元本構(gòu)關系方法將其等效為線性子結(jié)構(gòu)，通過所有子結(jié)構(gòu)合并計算求解整體非線性動力學方程，用較小計算成本獲得非線性體系既滿足工程應用且較精確的動力解。

2結(jié)構(gòu)剛度矩陣提取方法

用商業(yè)軟件進行各種結(jié)構(gòu)動、靜力分析可據(jù)不同目的提取剛度、質(zhì)量、阻尼及荷載矩陣，以便利用各類矩陣做后期分析與處理。對線彈性結(jié)構(gòu)可利用編程提取，方法較簡單；但對非線性彈塑性結(jié)構(gòu)各類矩陣提取較困難，大型商業(yè)軟件一般無直接提取方法。致采用變化的結(jié)構(gòu)特性矩陣研究結(jié)構(gòu)某些階段乃至整個階段受力特點造成困難。本文利用ANSYS二次開發(fā)工具嘗試進行特性矩陣提取編程，其中剛度矩陣提取為關注重點。

在ANSYS有限元程序中提取整體剛度矩陣方法主要是HBMAT命令法及超單元(Super-element)法。HBMAT命令為ANSYS提取整體剛度矩陣的直接內(nèi)部提取方法，但其僅適用于線彈性分析，不適用非線性分析，且該命令采用索引存儲方法的稀疏矩陣并以Harwell-Boeing格式存儲剛度矩陣下三角非零數(shù)據(jù)，并需后期借用MATLAB等程序編寫命令處理還原其矩陣形式。同時HBMAT命令法也有嚴重不足，即形成的剛度矩陣無序，對不同結(jié)構(gòu)單元網(wǎng)格較推算出矩陣某行某列的值與哪個節(jié)點自由度相關，且不能說明整體剛度如何從單元剛度矩陣對號入座組裝的，而此過程恰為關注內(nèi)容。超單元法同樣僅可提取線彈性整體剛度矩陣，基本步驟為，創(chuàng)建有限元模型并施加約束、定義分析類型為子結(jié)構(gòu)、定義輸入何種矩陣、選擇并定義所有節(jié)點為主自由度、求解并列出矩陣。所列整體剛度矩陣為全部元素(全矩陣)，按行順序分別列出各列元素數(shù)值，但結(jié)構(gòu)節(jié)點較多時，其數(shù)據(jù)量龐大，且后期提取數(shù)據(jù)需大量人工處理，非常繁瑣導致不方便使用。

由于HBMAT命令法及超單元法在提取整體剛度矩陣的應用上均有不便之處，故本文基于FORTRAN語言的UPFs二次開發(fā)工具在ANSYS有限元程序中提取矩陣。其核心思想為通過編寫外部用戶程序從ANSYS子空間計算方法模態(tài)分析結(jié)果File.Full的二進制文件中提取特性矩陣各行、列非零值，編寫MATLAB程序按節(jié)點編號由小到大順序?qū)μ柸胱M裝形成完整矩陣形式。ANSYS二次開發(fā)環(huán)境為Compaq Fortran 6.5，其中主文件為Matrix-extraction. For，其余Matrix-output.F90用于矩陣輸出，Binlib.Lin為ANSYS提供庫文件，Binlib.Dll為動態(tài)鏈接庫文件。運行編譯后形成可執(zhí)行Matrix-output.Exe文件即獲得質(zhì)量矩陣(Mass.Matrix)、剛度矩陣(Stiffness.Matrix)文件。

3非線性等效剛度矩陣處理方法

考慮結(jié)構(gòu)在非線性塑性階段提取各時刻剛度矩陣，可近似認為該非線性塑性本構(gòu)關系由分段逐步遞減且在小段Δt時間內(nèi)線性本構(gòu)關系組合。對進入非線性階段結(jié)構(gòu)，在小段Δt時間內(nèi)提取t及t+Δt時刻各單元節(jié)點應力應變。以2D平面應變問題為例進行推導，每個單元每個節(jié)點平面應變的物理方程為

t時刻

(1)

t+Δt時刻

(2)

則Δt時段應變差為

(3)

式(3)為各向同性的物理方程。若為提取剛度矩陣修改E，G則變成各項異性，物理增量方程應為

(4)

式(4)可簡化為

(5)

(6)

ANSYS每時刻Step運行結(jié)束后調(diào)用MATLAB程序時須設置只有在該Step調(diào)用的MATLAB程序分析完成后才可繼續(xù)運行ANSYS程序，即ANSYS在調(diào)用MATLAB程序未計算完畢時，須等待而不能繼續(xù)進行運算。需兩者同時運行并建立Flag文件，通過兩者中讀其內(nèi)容判斷對方是否運行。兩者若運行完一個Step，改變Flag，告訴對方自己當前運行結(jié)束，對方可繼續(xù)運行，否則必須等待。

4局部非線性SSI體系的CMS應用

為檢驗局部非線性下約束模態(tài)綜合法及勢能判據(jù)截斷準則的有效性，擬對地基土-高層剪力墻結(jié)構(gòu)相互作用體系進行動力分析。計算模型見圖1。上部為10層剪力墻結(jié)構(gòu)，層高3.6 m，寬24 m；支撐于2層箱型基礎上。土體設為分層土，地質(zhì)參數(shù)見表1。土體區(qū)域沿承臺兩側(cè)寬度各取100 m，沿深度取30 m?？蚣芗袅εc箱型基礎阻尼比為0.05，土阻尼比設為0.2。結(jié)構(gòu)模型兩端自由端，底端固定端。對整體結(jié)構(gòu)進行Taft波、EL-centro波及天津人工波動力時程分析。篇幅所限，以天津人工波為例進行分析。圖2為ANSYS直接計算法中地震波作用下其下部土體塑性區(qū)變形云圖，表現(xiàn)出局部非線性特性。為比較約束模態(tài)綜合法的優(yōu)越性及精確度，將計算結(jié)果與ANSYS直接計算法進行對比分析。

圖1　地基土-高層剪力墻結(jié)構(gòu)體系計算模型 Fig.1 Calculation model of soil-structure interaction system

編號厚度/m彈模/Pa泊松比μ粘聚力/Pa摩擦角Φ/(°)濕容重/(kg·m-3)17.12.00E70.21425014.51940212.03.48E70.21613917.31915310.95.00E70.21850023.92036

圖2　地震波作用下整體模型塑性區(qū)變形云圖 Fig.2 Von misis equivalent plastic strain contour

圖3　子結(jié)構(gòu)示意圖 Fig.3 Substructure division

圖4　子結(jié)構(gòu)1截斷參數(shù)判定 Fig.4 Substructure 1 truncation parameter determination

本文采用逐步提取剛度矩陣法可方便地通過不同時刻特征值方程求得結(jié)構(gòu)頻率變化趨勢。圖8為天津人工波作用下整體結(jié)構(gòu)前5階頻率變化過程。由圖8看出，結(jié)構(gòu)在彈性階段頻率保持不變，進入塑性階段時呈階梯狀降低。

圖5　位移時程曲線對比 Fig.5 Displacement curve comparison

圖6 速度時程曲線對比Fig.6Velocitycurvecomparison圖7 加速度時程曲線對比Fig.7Accelerationcurvecomparison圖8 地震波作用下前3階頻率變化曲線Fig.8First3orderfrequencycurvesunderearthquakeaction

天津人工地震波動力非線性分析前3階頻率在整個過程的某些時間點變化見表2，可見當?shù)卣鸩ㄟ_到峰值時基頻降低程度最大，趨于平穩(wěn)時各階頻率保持恒定。

表2　各階頻率變化情況 (Hz)

5結(jié)論

(1)本文據(jù)土-結(jié)構(gòu)相互作用存在局部非線性特點用線性-非線性混合約束模態(tài)綜合法進行分析。結(jié)果表明，用該方法與有限元直接法所得響應曲線基本吻合，驗證了在彈塑性階段利用分段等效線性化手段及ANSYS二次開發(fā)程序文件提取矩陣的正確性與可行性。與有限元直接法相比，其計算成本大幅度降低，是求解土-結(jié)構(gòu)動力相互作用問題的有效方法。

(2)本文方法在分析過程中對局部非線性部分提取剛度、質(zhì)量、阻尼矩陣為采用變化的特性矩陣對結(jié)構(gòu)進入塑性階段深層次分析研究提供有效方法。

參考文獻

[1]王菲. 地基土-高層建筑相互作用的動態(tài)子結(jié)構(gòu)法[D].天津：天津大學,2010.

[2]白建方, 樓夢麟. 基于動力子結(jié)構(gòu)方法的場地地震反應分析方法[J]. 震災防御技術, 2008,3(2)：145-154.

BAI Jian-fang, LOU Meng-lin. The dynamic sub-structure method for seismic response of irregular topography[J]. Technology for Earth-quake Disaster Prevention, 2008, 3(2)：145-154.

[3]曲哲, 葉列平. 計算結(jié)構(gòu)非線性地震峰值響應的等價線性化模型[J]. 工程力學, 2011,28(10)：93-100.

QU Zhe, YE Lie-ping. An equivalent linear model to estimate maximum inelastic seismic responses of structural systems[J]. Engineering Mechanics,2011, 28(10) ：93-100.

[4]閆相橋, 杜善義, 王鐸. 材料非線性有限元分析中組織結(jié)構(gòu)剛度陣的子結(jié)構(gòu)法[J]. 哈爾濱工業(yè)大學學報,1989,12(6)： 100-103.

YAN Xiang-qiao, DU Shan-yi, WANG Duo. An effective method of assembling the structural stiffness matrix in material nonlinear finite element analyses[J]. Journal of Harbin Institute of Technology,1989,12(6):100-103.

[5]柳國環(huán), 陸新征. 基巖地震譜與地震動位移輸入的土-結(jié)構(gòu)相互作用(SSI)計算模型改進[J]. 巖石力學與工程學報, 2011,30(5):84-89.

LIU Guo-huan, LU Xin-zheng. Spectra of bedrock earthquake motion and improvement of soil-structure interaction(SSI) calculation model for seismic displacement inputting[J]. Chinese Journal of Rock Mechanics and Engineering,2011,30(5):84-89.

[6]應祖光, 葉淑琴, 金林. 基于固定界面子結(jié)構(gòu)模態(tài)的頻響函數(shù)精確綜合法[J]. 振動與沖擊,2010,29(3):132-133.

YING Zu-guang, YE Shu-qin, JIN Lin. Exact frequency response function synthesis method using interface fixed substructure modes[J]. Journal of Vibration and Shock, 2010, 29(3):132-133.