土-結構非線性相互作用混合約束模態實施方法

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

郵箱：jiangxinliang@126.com

土-結構非線性相互作用混合約束模態實施方法

姜忻良，張海順

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

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

關鍵詞：土-結構相互作用；約束模態綜合法；勢能判據截斷準則；局部非線性；彈塑性剛度矩陣

收稿日期：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

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

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

1基本理論及方法闡述

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

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

2結構剛度矩陣提取方法

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

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

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

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

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

t時刻

(1)

t+Δt時刻

(2)

則Δt時段應變差為

(3)

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

(4)

式(4)可簡化為

(5)

(6)

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

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

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

圖1　地基土-高層剪力墻結構體系計算模型 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　地震波作用下整體模型塑性區變形云圖 Fig.2 Von misis equivalent plastic strain contour

圖3　子結構示意圖 Fig.3 Substructure division

圖4　子結構1截斷參數判定 Fig.4 Substructure 1 truncation parameter determination

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

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

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

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

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

5結論

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

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

參考文獻

[1]王菲. 地基土-高層建筑相互作用的動態子結構法[D].天津：天津大學,2010.

[2]白建方, 樓夢麟. 基于動力子結構方法的場地地震反應分析方法[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]曲哲, 葉列平. 計算結構非線性地震峰值響應的等價線性化模型[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]閆相橋, 杜善義, 王鐸. 材料非線性有限元分析中組織結構剛度陣的子結構法[J]. 哈爾濱工業大學學報,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]柳國環, 陸新征. 基巖地震譜與地震動位移輸入的土-結構相互作用(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]應祖光, 葉淑琴, 金林. 基于固定界面子結構模態的頻響函數精確綜合法[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.