姜忻良，張海順

(1. 天津大學(xué)建筑工程學(xué)院，天津 300072；2. 濱海土木工程結(jié)構(gòu)與安全教育部重點(diǎn)實(shí)驗(yàn)室(天津大學(xué))，天津 300072)

動(dòng)力時(shí)程分析中彈塑性剛度矩陣的提取方法

姜忻良1,2，張海順1,2

(1. 天津大學(xué)建筑工程學(xué)院，天津 300072；2. 濱海土木工程結(jié)構(gòu)與安全教育部重點(diǎn)實(shí)驗(yàn)室(天津大學(xué))，天津 300072)

利用ANSYS 二次開發(fā)工具UPFs編寫成可執(zhí)行程序文件，通過與MATLAB進(jìn)行聯(lián)立計(jì)算，對(duì)線彈性和非線性彈塑性本構(gòu)模型的結(jié)構(gòu)分別實(shí)現(xiàn)了特性矩陣的近似提?。畬?duì)線性本構(gòu)關(guān)系模型的各特性矩陣提取方法進(jìn)行闡述推導(dǎo)和驗(yàn)證，求得的各階頻率和動(dòng)態(tài)響應(yīng)時(shí)程曲線與 ANSYS直接計(jì)算法吻合．對(duì)于非線性彈塑性本構(gòu)模型，將非線性本構(gòu)關(guān)系在全局過程中的小段時(shí)間內(nèi)進(jìn)行等效線性化處理后提取等效剛度矩陣，進(jìn)行靜動(dòng)力分析后所得各種響應(yīng)時(shí)程曲線與 ANSYS直接計(jì)算法所得曲線基本吻合．結(jié)果表明了分段等效線性化模型的合理性和ANSYS二次開發(fā)程序提取等效特性矩陣的可行性．

ANSYS二次開發(fā)；彈塑性剛度矩陣；分段線性化；等效彈性模量；稀疏矩陣；模態(tài)分析；動(dòng)力時(shí)程分析

在采用商業(yè)軟件進(jìn)行各種結(jié)構(gòu)動(dòng)、靜力分析時(shí)，研究人員根據(jù)不同目的而需要對(duì)剛度、質(zhì)量、阻尼和荷載矩陣進(jìn)行提取，以方便利用上述各類矩陣做后期分析與處理[1-3]．其中對(duì)于線彈性結(jié)構(gòu)，研究人員可利用編程手段進(jìn)行提取，方法較為簡單[4-5]；但對(duì)于非線性彈塑性結(jié)構(gòu)，上述各類矩陣的提取較為困難，目前大型商業(yè)軟件一般無直接簡便的提取方法．這給采用變化的結(jié)構(gòu)特性矩陣來研究結(jié)構(gòu)某些階段或整個(gè)階段的受力特點(diǎn)帶來困難．

本文利用ANSYS二次開發(fā)工具手段來嘗試進(jìn)行提取編程，而其中剛度矩陣(包括單元?jiǎng)偠染仃嚒⒃冀Y(jié)構(gòu)剛度矩陣和結(jié)構(gòu)剛度矩陣等)的提取是最為研究人員所關(guān)心的，本文主要對(duì)提取結(jié)構(gòu)剛度矩陣進(jìn)行嘗試分析．對(duì)線性本構(gòu)關(guān)系模型的剛度矩陣基于ANSYS二次開發(fā)提取方法進(jìn)行闡述推導(dǎo)和驗(yàn)證，分別由特征值方程計(jì)算求得頻率和由地震波動(dòng)力時(shí)程分析求得響應(yīng)曲線，并與 ANSYS直接計(jì)算法所得到的頻率和響應(yīng)曲線進(jìn)行對(duì)比分析，以驗(yàn)證所二次開發(fā)的程序文件提取矩陣的可行性．

1 ANSYS二次開發(fā)工具

用于ANSYS二次開發(fā)的工具主要有4個(gè)，即APDL、UPFs(user programmable features)、UIDL和Tcl/Tk，使用以上工具可以建立新的材料模型，構(gòu)建新的單元類型，參數(shù)化建模，優(yōu)化分析，構(gòu)建流程化的 ANSYS分析平臺(tái)，建立符合用戶專業(yè)需求的ANSYS用戶界面等[6-8]．

本文主要依據(jù)APDL和UPFs兩種工具，其中APDL應(yīng)用比較普遍，不再贅述介紹，而UPFs是在ANSYS提供的 Fortran源代碼的基礎(chǔ)上，修改其用戶可編程子程序和函數(shù)，從源代碼層次上對(duì)ANSYS進(jìn)行二次開發(fā)的工具．用戶需要在響應(yīng)的Fortran語言編譯器的支持下，將編譯修改后的源代碼與ANSYS庫相連而形成用戶所需的ANSYS可執(zhí)行文件來進(jìn)行后續(xù)操作．

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

在ANSYS有限元程序中提取整體剛度矩陣方法主要是 HBMAT命令法和超單元(super-element)法．HBMAT命令是ANSYS提取整體剛度矩陣的直接內(nèi)部提取方法，但該命令僅適用于線彈性分析，對(duì)非線性分析不適用，且該命令采用索引存儲(chǔ)的稀疏矩陣并以Harwell-Boeing格式來存儲(chǔ)剛度矩陣下三角的非零數(shù)據(jù)，并需要后期借用 MATLAB等程序編寫命令處理來還原其矩陣形式．但是同時(shí)HBMAT命令法也有嚴(yán)重不足，首先形成的剛度矩陣是無序的，對(duì)于不同結(jié)構(gòu)的單元網(wǎng)格，研究人員很難推算出矩陣的某一行某一列的值是與哪個(gè)節(jié)點(diǎn)自由度所相關(guān)的，且不能說明整體剛度是如何從單元?jiǎng)偠染仃嚒皩?duì)號(hào)入座”組裝來的，而此過程卻恰恰是研究人員所關(guān)心的內(nèi)容．

超單元法同樣僅可提取線彈性整體剛度矩陣，基本步驟是：創(chuàng)建有限元模型并施加約束，定義分析類型為子結(jié)構(gòu)，定義輸入何種矩陣，選擇并定義所有節(jié)點(diǎn)為主自由度，求解并列出矩陣．所列出的整體剛度矩陣為全部元素(全矩陣)，按行的順序分別列出各列元素?cái)?shù)值，但其問題是當(dāng)結(jié)構(gòu)節(jié)點(diǎn)較多時(shí)，數(shù)據(jù)量非常龐大，且后期提取數(shù)據(jù)需要大量人工手動(dòng)處理，非常繁瑣導(dǎo)致使用不方便．

上述HBMAT命令法和超單元法在提取整體剛度矩陣的應(yīng)用上均有不便之處，所以本文嘗試基于Fortran語言的UPFs工具在ANSYS有限元程序中進(jìn)行二次開發(fā)來提取矩陣．其核心思想就是通過編寫外部用戶程序直接從ANSYS子空間計(jì)算方法的模態(tài)分析結(jié)果的File.Full二進(jìn)制文件中提取矩陣的各行各列非零值，然后編寫 MATLAB程序使其有序地按照節(jié)點(diǎn)編號(hào)由小到大的順序“對(duì)號(hào)入座”組裝成為完整的矩陣形式．

ANSYS二次開發(fā)環(huán)境為Compaq Fortran 6.5，其中主文件為 Matrix-extraction.For，其余 Matrixoutput.F90用于矩陣輸出，Binlib.Lin為 ANSYS提供庫文件，Binlib.Dll為動(dòng)態(tài)鏈接庫文件．運(yùn)行編譯后而形成Matrix-output.Exe文件即可得到質(zhì)量矩陣(mass matrix)和剛度矩陣(stiffness matrix)文件．整個(gè)分析的具體的流程如下．

(1) 線彈性本構(gòu)模型ANSYS-MATLAB聯(lián)立計(jì)算．對(duì)于線彈性本構(gòu)模型的結(jié)構(gòu)，首先提取模型的各個(gè)節(jié)點(diǎn)編號(hào)，進(jìn)行ANSYS模態(tài)分析后按節(jié)點(diǎn)編號(hào)順序在MATLAB中對(duì)號(hào)入座生成初始剛度矩陣和質(zhì)量矩陣，并由此計(jì)算Rayleigh阻尼矩陣和荷載矩陣．

(2) 彈塑性本構(gòu)模型ANSYS-MATLAB聯(lián)立計(jì)算．對(duì)于非線性彈塑性本構(gòu)模型的結(jié)構(gòu)，每個(gè)時(shí)刻計(jì)算完畢后判斷結(jié)構(gòu)是否進(jìn)入塑性階段，若未進(jìn)入塑性階段則仍利用方法(1)的線彈性計(jì)算步驟；若已經(jīng)進(jìn)入塑性階段，則首先提取此時(shí)的各個(gè)塑性單元的應(yīng)力和應(yīng)變值，通過由公式推導(dǎo)求出的等效彈性模量和等效剪切模量，并賦予該單元新的材料屬性；重新進(jìn)行ANSYS模態(tài)分析后，通過MATLAB程序提取其等效剛度矩陣、質(zhì)量矩陣、阻尼矩陣和荷載矩陣；上述步驟均完成之后再進(jìn)行動(dòng)力時(shí)程分析從而結(jié)束一個(gè)時(shí)刻的計(jì)算．隨后將其本構(gòu)模型恢復(fù)至初始狀態(tài)，利用 ANSYS的重啟命令繼續(xù)下一時(shí)刻的計(jì)算，如此往復(fù)進(jìn)行．

3 線性模型的特性矩陣提取方法的驗(yàn)證

為檢驗(yàn)上述ANSYS二次開發(fā)方法所提取的剛度矩陣和質(zhì)量矩陣的正確性，有必要做一般性驗(yàn)證分析．對(duì)一簡單彈性平面應(yīng)變模型(見圖 1)進(jìn)行模態(tài)分析和天津人工地震波(見圖 2)作用下的動(dòng)力分析．

圖1 網(wǎng)格模型Fig.1 Unit grid model

圖2 天津人工地震波Fig.2 Tianjin artificial wave

模型采用平面應(yīng)變 PLANE42單元，彈性模量E＝80,MPa，泊松比μ＝0.3，密度ρ＝1,750,kg/m3，尺寸為 5,m×12,m，底部固定．上述分析均采用ANSYS直接計(jì)算法和ANSYS-MATLAB聯(lián)立法計(jì)算并進(jìn)行了對(duì)比分析．對(duì)比分析內(nèi)容包括頻率、動(dòng)態(tài)響應(yīng)時(shí)程曲線和計(jì)算時(shí)間．

用ANSYS直接進(jìn)行模態(tài)分析計(jì)算并提取所有模態(tài)頻率，再利用 ANSYS二次開發(fā)程序所提取的剛度矩陣K和質(zhì)量矩陣M在MATLAB中進(jìn)行特征值方程運(yùn)算，求得各階頻率進(jìn)行對(duì)比，見表 1．從表 1中可以看出兩種方法計(jì)算所得的各階頻率相同，并且采用二次開發(fā)程序所運(yùn)行的計(jì)算時(shí)間大為減少，表明此 ANSYS二次開發(fā)提取矩陣運(yùn)算的正確性和高效性．

動(dòng)力時(shí)程分析中對(duì)整體結(jié)構(gòu)施加天津人工地震波，提取其頂端中點(diǎn)的位移曲線待用，此過程由ANSYS直接計(jì)算；然后，將 ANSYS和 MATLAB兩個(gè)程序聯(lián)立，同樣利用二次開發(fā)方法提取特性矩陣(剛度、質(zhì)量、阻尼和荷載矩陣)，求得其頂端中點(diǎn)位移時(shí)程曲線并與 ANSYS的計(jì)算結(jié)果對(duì)比，如圖3所示．由圖 3可以看出，提取的特性矩陣所計(jì)算的頂端中點(diǎn)位移時(shí)程曲線與ANSYS直接計(jì)算的結(jié)果吻合；速度、加速度曲線為位移曲線對(duì)時(shí)間求導(dǎo)所得，故也吻合．此外，二次開發(fā)程序所運(yùn)行的動(dòng)力時(shí)程分析的計(jì)算時(shí)間也比ANSYS直接計(jì)算法的時(shí)間大幅度縮減，這也同樣說明了上述矩陣提取方法的正確性和高效性．

表1 兩種方法計(jì)算所得模態(tài)頻率Tab.1 Modal frequency of two methods

圖3 天津人工地震波矩陣驗(yàn)證Fig.3 Verification of matrix of Tianjin artificial wave

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

第3節(jié)中的分析驗(yàn)證了ANSYS二次開發(fā)在線彈性模型下所提取的矩陣真實(shí)有效．若考慮結(jié)構(gòu)在非線性塑性階段提取其各個(gè)時(shí)刻的剛度矩陣，可以近似地認(rèn)為非線性塑性本構(gòu)關(guān)系是分段逐步遞減的在小時(shí)間段 Δt內(nèi)的線性本構(gòu)關(guān)系的組合．因此，若要在動(dòng)力時(shí)程分析的塑性階段生成剛度矩陣，就需要將 Δt時(shí)間內(nèi)視為線性本構(gòu)關(guān)系．

對(duì)于進(jìn)入非線性階段的結(jié)構(gòu)，在小時(shí)間段 Δt內(nèi)，提取t和t+ Δt時(shí)刻的各個(gè)單元節(jié)點(diǎn)的x、y方向的正應(yīng)力 σx、σy和剪應(yīng)力τxy，以及與之對(duì)應(yīng)的 x、y方向的正應(yīng)變 εx、εy和剪應(yīng)變 γxy．對(duì)于平面應(yīng)變問題，依據(jù)如下公式進(jìn)行推導(dǎo)．各單元的每一個(gè)節(jié)點(diǎn)的平面應(yīng)變的物理方程為

t時(shí)刻

t＋ Δt時(shí)刻

則 Δt時(shí)段的應(yīng)變差為

求解式(5)得到該單元節(jié)點(diǎn)的等效彈性模量和等效剪切模量，即

式(3)為各向同性的物理方程，若為了提取剛度矩陣而修改E、G則變成各項(xiàng)異性，則物理增量方程應(yīng)為

式(4)可以簡化為

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

5 非線性模型特性矩陣提取方法驗(yàn)證

第 4節(jié)推導(dǎo)出非線性塑性階段分段等效剛度矩陣的計(jì)算提取方法，下面分別應(yīng)用Pushover靜力非線性分析和天津人工地震波的動(dòng)力非線性分析方法來進(jìn)行該方法的驗(yàn)證．

模型同樣采用圖 1所示的結(jié)構(gòu)，其非線性本構(gòu)關(guān)系采用理想彈塑性的 DP模型，黏聚力 C為20,kPa，內(nèi)摩擦角和膨脹角均為 30°．在 Pushover靜力非線性分析中，假定在某一荷載步時(shí)，某一單元的最高等效塑性應(yīng)變?cè)隽砍^該單元前一步等效塑性應(yīng)變峰值的 15%時(shí)，認(rèn)為此結(jié)構(gòu)進(jìn)入強(qiáng)非線性階段，接近破壞而停止加載[9-10]．模型的 von Misis等效塑性應(yīng)變?cè)茍D如圖4所示．

將ANSYS和MATLAB兩個(gè)程序?qū)樱诿恳粋€(gè)荷載步結(jié)束后進(jìn)行后處理分析．依據(jù)式(6)計(jì)算每個(gè)單元的等效彈性模量和等效剪切模量隨后賦予該單元材料屬性，繼而進(jìn)行模態(tài)分析并對(duì)整體結(jié)構(gòu)提取其剛度矩陣與質(zhì)量矩陣，再調(diào)用MATLAB程序通過特征值方程依次求出前3階頻率．全過程的前 3階頻率變化趨勢(shì)如圖 5所示，可以看出結(jié)構(gòu)在進(jìn)入塑性階段后，各階頻率逐漸降低，說明其剛度在逐漸減?。?/p>

圖4 模型的von Misis等效塑性應(yīng)變?cè)茍D(Pushover)Fig.4 von Misis equivalent plastic strain contour of the model(Pushover)

圖5 前3階頻率的變化曲線(Pushover)Fig.5 Variation curves of the 1st three order frequencies (Pushover)

依據(jù)上述各個(gè)時(shí)刻的矩陣，利用 MATLAB通過 Hooke定律 Ku =F，求得其頂端中點(diǎn)位移時(shí)程曲線和底部最大塑性區(qū)節(jié)點(diǎn)位移曲線，與 ANSYS直接計(jì)算法的曲線進(jìn)行對(duì)比，結(jié)果如圖6所示．

圖6 模型的位移曲線(Pushover)Fig.6 Displacement curves of the model(Pushover)

可以看出圖 6中兩種方法的位移曲線比較吻合，說明在非線性彈塑性階段所提取的等效彈性模量和等效剪切模量較為準(zhǔn)確，繼而說明式(6)推導(dǎo)的在 Δt時(shí)間內(nèi)非線性彈塑性本構(gòu)模型可以轉(zhuǎn)變?yōu)榉侄蔚刃Ь€性化模型進(jìn)行處理，也驗(yàn)證了本文ANSYS二次開發(fā)程序的可行性．

圖 7中為天津人工地震波動(dòng)力分析的結(jié)果．由圖 7可見在地震波作用下，結(jié)構(gòu)底部兩端均進(jìn)入塑性階段，與實(shí)際情況較符合．此外，利用 ANSYS直接計(jì)算法和ANSYS-MATLAB聯(lián)立法得到的結(jié)構(gòu)頂部中點(diǎn)位移、速度、加速度時(shí)程曲線見圖 8．這些動(dòng)力響應(yīng)(位移、速度、加速度)時(shí)程曲線基本上吻合，說明采用等效彈性模量和等效剪切模量所提取的結(jié)構(gòu)逐漸變化的剛度矩陣較為準(zhǔn)確.

圖7 模型的von Misis等效塑性區(qū)應(yīng)變?cè)茍D(地震波)Fig.7 von Misis equivalent plastic strain contour (earthquake wave)

圖8 兩種方法的頂部中點(diǎn)動(dòng)力響應(yīng)時(shí)程曲線對(duì)比Fig.8 Contrast of dynamic response time history curves of top midpoint between two methods

圖 9為地震波作用下結(jié)構(gòu)前 3階頻率的變化過程，可以看出結(jié)構(gòu)在彈性階段內(nèi)頻率保持不變，在進(jìn)入塑性階段時(shí)，結(jié)構(gòu)的頻率呈階梯狀降低．

圖9 地震波作用下前3階頻率的曲線Fig.9 Curves of the 1st three order frequencies under earthquake action

表 2和表 3匯集了上述 Pushover靜力非線性分析和天津人工地震波動(dòng)力非線性分析的前 3階頻率在整個(gè)過程的某些時(shí)間點(diǎn)的變化情況．從Pushover靜力非線性分析中可以看出，隨著荷載的逐步增加，結(jié)構(gòu)的前 3階頻率均逐步降低，而第 1階頻率降低最為突出．在天津人工波動(dòng)力非線性分析中，在地震波達(dá)到峰值的時(shí)刻前 3階頻率的降低速度較大，當(dāng)?shù)卣鸩ㄚ呌谄椒€(wěn)的時(shí)候其前 3階頻率保持恒定，同樣也是第1階頻率降低得最多．

表2 前3階頻率變化情況(Pushover)Tab.2 Variation of the 1st three order frequencies(Pushover)

表3 前3階頻率變化情況(地震波)Tab.3 Variation of the 1st three order frequencies(Earthquake wave)

6 結(jié) 語

采用分段等效線性化的處理方法，對(duì)非線性彈塑性本構(gòu)模型在 Pushover靜力方法與時(shí)程分析中進(jìn)行各種矩陣的提?。Y(jié)果表明在非線性彈塑性階段所提取并采用的等效彈性模量和等效剪切模量較為準(zhǔn)確，說明在 Δt時(shí)間段內(nèi)將非線性彈塑性本構(gòu)模型轉(zhuǎn)變?yōu)榉侄蔚刃Ь€性化模型的合理性，同時(shí)也驗(yàn)證了本文 ANSYS二次開發(fā)程序的可行性．

本文方法為采用變化的特性矩陣對(duì)結(jié)構(gòu)進(jìn)入塑性階段后深層次的分析研究提供了有效的方法．

［1］ 張令心，石 磊. 土-結(jié)構(gòu)相互作用地震反應(yīng)分析軟件及其二次開發(fā)[J]. 地震工程與工程振動(dòng)，2006，26(3)：225-227.

Zhang Lingxin，Shi Lei. Software of seismic soil structure interaction analysis and its re-developing[J]. Earthquake Engineering and Engineering Vibration，2006，26(3)：225-227(in Chinese).

［2］ 劉艷萍，楊新華，楊文兵. 預(yù)應(yīng)力鋼筋混凝土局部有限元分析的 ANSYS二次開發(fā)[J]. 華中科技大學(xué)學(xué)報(bào)(城市科學(xué)版)，2005，22(增)：87-90.

Liu Yanping，Yang Xinhua，Yang Wenbing. ANSYS secondary development for the finite element analysis of prestressed reinforced concrete structures[J]. Journal of Huazhong University of Science and Technology(Urban Science Edition)，2005，22(Suppl)：87-90(in Chinese).

［3］ 王一功，楊佑發(fā). 針對(duì)場(chǎng)地地震反應(yīng)分析的 ANSYS二次開發(fā)[J]. 地震工程與工程振動(dòng)，2004，24(2)：42-45.

Wang Yigong，Yang Youfa. A secondary development of ANSYS for site earthquake response[J]. Earthquake Engineering and Engineering Vibration，2004，24(2)：42-45(in Chinese).

［4］ 劉光棟，王解君，何放龍. 空間梁單元的幾何非線性剛度矩陣的分解形式[J]. 湖南大學(xué)學(xué)報(bào)，1992，19(1)：60-71.

Liu Guangdong，Wang Jiejun，He Fanglong. Resolvedformulation of geometrical nonlinear stiffness matrix for three-dimensional beam element[J]. Journal of Hunan University，1992，19(1)：60-71(in Chinese).

［5］ 劉齊茂，燕柳斌. 基于 Newmark 法敏度計(jì)算的剛架結(jié)構(gòu)動(dòng)力優(yōu)化[J]. 工程力學(xué)，2010，27(3)：145-154.

Liu Qimao，Yan Liubin. Dynamic optimization of frame structures using sensitivity calculation based on newmark method[J]. Engineering Mechanics，2010，27(3)：145-154(in Chinese).

［6］ 李 妍，吳 斌，歐進(jìn)萍. 彈塑性結(jié)構(gòu)等效線性化方法的對(duì)比研究[J]. 工程抗震與加固改造，2005，27(1)：1-6.

Li Yan，Wu Bin，Ou Jinping. Comparison of equivalent linearization methods for inelastic structures[J]. Earthquake Resistant Engineering and Retrofitting，2005，27(1)：1-6(in Chinese).

［7］ 曲 哲，葉列平. 建筑結(jié)構(gòu)彈塑性地震響應(yīng)計(jì)算的等價(jià)線性化法研究[J]. 建筑結(jié)構(gòu)學(xué)報(bào)，2010，31(9)：95-102.

Qu Zhe，Ye Lieping. Equivalent linear analysis in estimating nonlinear seismic responses of building structures[J]. Journal of Building Structures，2010，31(9)：95-102(in Chinese).

［8］ 程光煜，葉列平. 基于等效線性化方法彈塑性

SDOF系統(tǒng)能量譜的研究[J]. 建筑結(jié)構(gòu)，2007，37(8)：74-77. Cheng Guangyu，Ye Lieping. Estimation of input energy spectra of inelastic SDOF systems with equivalent linear method[J]. Building Structures，2007，37(8)：74-77(in Chinese).

［9］ 白建方. 復(fù)雜場(chǎng)地土層地震反應(yīng)分析的并行有限元方法[D]. 上海：同濟(jì)大學(xué)土木工程學(xué)院，2007.

Bai Jianfang. Parallel Finite Element Method for Seismic Response Analysis of Irregular Site[D]. Shanghai：School of Civil Engineering，Tongji University，2007(in Chinese).

［10］ 張國棟. 土與結(jié)構(gòu)相互作用體系隨機(jī)地震反應(yīng)分析[D]. 武漢：武漢大學(xué)土木建筑工程學(xué)院，2004.

Zhang Guodong. Stochastic Seismic Response Analysis for Soil-Structure Interaction System[D]. Wuhan：School of Civil Engineering，Wuhan University，2004(in Chinese).

(責(zé)任編輯：樊素英)

Extraction Method of Elastic-Plastic Stiffness Matrix in Dynamic Time History Analysis

Jiang Xinliang1,2，Zhang Haishun1,2
(1. School of Civil Engineering，Tianjin University，Tianjin 300072，China；2. Key Laboratory of Coastal Civil Engineering Structure and Safety of Ministry of Education (Tianjin University)，Tianjin 300072，China)

Executable program for feature matrices extraction was written by using ANSYS secondary development tools UPFs. With MATLAB simultaneous calculation, the feature matrices of the linear elastic and nonlinear elastic-plastic constitutive model were respectively extracted. For the linear constitutive model, the extraction method of those feature matrices were derived and affirmed. All order frequencies and response curves were coincided with ANSYS direct calculation method on the linear constitutive model. Within short time of the entire period，the nonlinear constitutive model was transformed through equivalent linearization for extracting the equivalent stiffness matrix. The static and dynamic response time history curves were basically coincided with the ANSYS direct calculation method. Results show that the piecewise equivalent linearization is rational and the ANSYS secondary development program is feasible and accurate.

ANSYS secondary development；elastic-plastic stiffness matrix；piecewise linearization；equivalent elastic modulus；sparse matrix；modal analysis；dynamic time history analysis

TU317.1

A

0493-2137(2015)04-0355-07

10.11784/tdxbz201309007

2013-09-02；

2014-03-13．

國家自然科學(xué)基金資助項(xiàng)目(51178308，51278335)．

姜忻良（1951— ），男，博士，教授，jiangxinliang@126.com．

張海順，zhangyibiao_0216@163.com．