靜動力彈塑性分析的桿系離散單元計算理論研究

許玲玲，葉繼紅

(1. 河南工業(yè)大學(xué)土木工程學(xué)院，鄭州 450001；2. 中國礦業(yè)大學(xué)，江蘇省土木工程環(huán)境災(zāi)變與結(jié)構(gòu)可靠性重點實驗室，徐州 221116；3. 中國礦業(yè)大學(xué)，徐州市工程結(jié)構(gòu)火災(zāi)安全重點實驗室，徐州 221116)

以網(wǎng)殼結(jié)構(gòu)為代表的大跨空間結(jié)構(gòu)的靜動力彈塑性分析一直是研究熱點[1-3]。對于多由梁、桿等長細構(gòu)件組成的該類結(jié)構(gòu)，有限單元法的彈塑性分析方法有塑性鉸法、精細塑性鉸法、塑性區(qū)法等[4-5]。前兩種是以內(nèi)力為基本未知量建立屈服函數(shù)，后者則是以應(yīng)力作為判斷彈塑性狀態(tài)的依據(jù)。塑性鉸法[6]作為一種簡化算法，無法模擬截面上及沿桿長的塑性分布和發(fā)展，計算精度不高。塑性區(qū)法[7]將桿件劃分為多個區(qū)段，在區(qū)段間的截面上考慮塑性沿截面和桿長的發(fā)展過程，一般被認為是精確算法，然而，因其計算成本高而較少應(yīng)用于大型復(fù)雜結(jié)構(gòu)分析。精細塑性鉸法兼顧了上述兩種方法的優(yōu)點，即在計算成本允許的條件下，近似地考慮塑性對單元剛度的削弱。Liew 和Chen 等[8]通過穩(wěn)定插值函數(shù)模擬結(jié)構(gòu)的幾何非線性行為，在此基礎(chǔ)上建立了精細塑性鉸模型考慮殘余應(yīng)力對單元剛度的影響；Kim 等[9]提出了基于位移有限元法的精細塑性鉸模型，并對空間鋼架的彈塑性行為進行了分析。然而，有限單元法處理彈塑性問題時，結(jié)構(gòu)進入塑性后，單元剛度矩陣會成為非常量，此時平衡條件需在變形后的構(gòu)形上建立，故需迭代求解非線性方程組[10]。

傳統(tǒng)顆粒離散單元[11]已發(fā)展40 多年，主要用于解決散體的彈性力學(xué)問題，其在結(jié)構(gòu)工程的研究成果集中在地震作用、沖擊荷載下筋混凝土或砌體結(jié)構(gòu)的破壞過程[12-13]?；诖耍~繼紅教授課題組[14-17]提出了適用于空間桿件結(jié)構(gòu)的桿系離散單元法，并率先進行了系統(tǒng)研究。研究成果表明桿系離散單元法在結(jié)構(gòu)的彈塑性行為模擬中具有獨特優(yōu)勢，即幾何非線性問題和動力響應(yīng)的求解自動包含在顆粒的運動控制方程中，是一個自然過程；桿系離散單元法中結(jié)構(gòu)的彈塑性行為僅與接觸點的內(nèi)力求解相關(guān)，無需改變顆粒的運動控制方程，基本的計算框架不會改變。齊念[14]首先建立了桿系離散單元塑性鉸法，該法假設(shè)塑性集中在接觸點處，一旦接觸點的內(nèi)力狀態(tài)達到極限屈服面即認為接觸截面完全進入塑性，其無法模擬接觸截面上及沿接觸單元(即等效梁)長的塑性分布和發(fā)展。為此，進一步建立了桿系離散單元纖維模型，將顆粒間接觸本構(gòu)模型的單根彈簧等效為分布彈簧，并將分布彈簧看作是接觸截面的若干根纖維，單根纖維的受力狀態(tài)用單軸的應(yīng)力-應(yīng)變關(guān)系進行描述。葉繼紅和張梅[16]在其基礎(chǔ)上進行了完善，提出了桿系離散單元塑性區(qū)法，并推導(dǎo)了接觸截面在三維應(yīng)力狀態(tài)下的彈塑性本構(gòu)模型。葉繼紅和王佳[17]基于CPU-GPU 異構(gòu)平臺，構(gòu)建了單元級并行、節(jié)點級并行的桿系離散單元計算框架，該并行算法加速比最高可達12.7 倍。

然而，分析上述文獻發(fā)現(xiàn)，對于桿系離散單元法的研究仍存在如下不足：1)現(xiàn)有文獻均假定切向彈簧僅用于描述純剪力引起的純剪切變形，這與桿件結(jié)構(gòu)常見的受力特征相違背，原因是桿件結(jié)構(gòu)通常長細比較大，可忽略剪切變形的影響，即根據(jù)彎曲梁理論可認為切向位移(即撓度)是由剪力產(chǎn)生的彎曲變形引起的，并不是由剪力產(chǎn)生的截面剪切變形引起的。故基于該假定推導(dǎo)出的切向接觸剛度系數(shù)無法用于桿件結(jié)構(gòu)問題的求解；2)桿系離散單元纖維模型和桿系離散單元塑性區(qū)法應(yīng)用于大型復(fù)雜結(jié)構(gòu)的動力彈塑性分析時較為耗時。

本文首先重新定義了桿系離散單元法中接觸本構(gòu)模型的切向彈簧，并嚴謹推導(dǎo)了面向梁單元的各方向上接觸單元剛度系數(shù)計算公式。其次，將精細塑性鉸法引入桿系離散單元法，借鑒其采用切線模量Et的概念近似考慮由殘余應(yīng)力引起的沿接觸單元(即等效梁)長的逐漸屈服效應(yīng)，接觸截面的逐漸屈服則由剛度退化函數(shù)μ表示，豐富了桿系離散單元法的彈塑性計算理論，并用一個數(shù)值算例驗證了所提彈塑性計算理論的優(yōu)勢和精確性。最后，采用桿系離散單元精細塑性鉸法對一個1/3.5 縮尺比的單層球面網(wǎng)殼振動臺試驗進行彈塑性分析，將其與試驗結(jié)果進行對比，進一步驗證本文方法的精度和穩(wěn)定性。

1 桿系離散單元計算理論的修正

任意空間網(wǎng)殼結(jié)構(gòu)的桿系離散單元模型如圖1所示。結(jié)構(gòu)被離散為一組有限的剛性顆粒，顆粒是定義結(jié)構(gòu)的位置坐標、質(zhì)量、邊界條件以及變形等的載體。顆粒間的連接即是接觸單元，根據(jù)接觸單元性質(zhì)建立的接觸力與接觸位移的關(guān)系即是顆粒間的接觸本構(gòu)模型。因此，桿系離散單元的基本理論主要包括顆粒運動方程的建立與求解、接觸本構(gòu)模型的構(gòu)建兩大部分，可見文獻[14 - 16]。所有顆粒的運動遵循牛頓第二定律，則任一顆粒的運動控制方程為：

圖1 空間網(wǎng)殼結(jié)構(gòu)離散模型Fig. 1 Discrete element model of a spatial reticulated shell structure

式中：M、I分別為整體坐標系下顆粒的質(zhì)量和轉(zhuǎn)動慣量；u¨ 、 θ¨分別為整體坐標系下顆粒的平動加速度和轉(zhuǎn)動加速度；F、M分別為整體坐標系下顆粒所受的力和力矩，包括內(nèi)力和外力。

顆粒所受內(nèi)力由與之發(fā)生接觸的顆粒構(gòu)成的接觸單元提供，然后，接觸單元處的內(nèi)力可通過力系平移定理轉(zhuǎn)移至該顆粒形心，而接觸點處的內(nèi)力增量則由接觸本構(gòu)模型計算[14]，其表達式為：

式中： Δf、 Δm為接觸點處內(nèi)力和內(nèi)力矩增量；K、Kθ分別為接觸本構(gòu)模型的平動和轉(zhuǎn)動剛度矩陣； Δu、 Δθ分別為局部坐標系下接觸點處的純平動位移和轉(zhuǎn)動位移，可通過剛體運動學(xué)求得，見文獻[14]。

圖2 為圖1 中任意兩相鄰顆粒A和B間的接觸本構(gòu)模型，該接觸本構(gòu)模型通過6 個彈簧表征顆粒間的接觸力和接觸力矩分量，每個彈簧對應(yīng)一個力或力矩類型，彈簧的變形即是接觸點C的位移和轉(zhuǎn)角增量。因此，如何合理而精確地推導(dǎo)各變量的接觸剛度系數(shù)即式(2)中的K和Kθ是保證桿系離散單元計算精度的關(guān)鍵所在。

圖2 面向梁單元的接觸本構(gòu)模型構(gòu)成Fig. 2 Beam element-oriented construction of contact constitutive model

現(xiàn)有文獻中均假定接觸本構(gòu)模型中的切向彈簧僅用于描述純剪力引起的純剪切變形，這與桿件結(jié)構(gòu)常見的受力特征相違背，故基于該假定推導(dǎo)出的切向接觸剛度系數(shù)無法用于桿件結(jié)構(gòu)問題的求解，而現(xiàn)有文獻給出的處理方式是直接假定一個數(shù)值，然后用大量算例進行驗證。簡而言之，已有成果中關(guān)于切向接觸剛度系數(shù)有如下困境，詳細推導(dǎo)得到的計算公式無法使用，而采用的切向接觸剛度系數(shù)又缺乏系統(tǒng)的推導(dǎo)過程。

本文以任一平面等截面直桿為例給出接觸剛度系數(shù)的推導(dǎo)過程。該直桿的截面構(gòu)造和材料特性為：截面積A、慣性矩I、彈性模量E、剪切模量G、截面切應(yīng)力不均勻系數(shù)k。將直桿離散為單排剛性顆粒，取模型中任意兩相鄰顆粒A和顆粒B，且把顆粒A和B間的接觸單元看作一根等效梁，該等效梁的長度L為兩顆粒形心的距離，截面構(gòu)造和材料特性與直桿一致。已有文獻[14]通過設(shè)置3 根獨立的彈簧表征相鄰顆粒間的軸力、剪力和彎矩，對應(yīng)的彈簧類型分別為法向彈簧、切向彈簧和轉(zhuǎn)動彈簧，如圖3 所示。其中，法向彈簧和轉(zhuǎn)動彈簧的剛度系數(shù)由公式推導(dǎo)獲得，可精確地描述顆粒間單向受拉壓變形(圖3(a))和由相對轉(zhuǎn)動引起的純彎曲行為(圖3(b))。切向彈簧描述的是由剪力引起顆粒間的純剪切變形(圖3(c))，然而，桿件結(jié)構(gòu)的長細比較大，通?？珊雎越孛婕羟凶冃蔚挠绊懀凑J為切向位移是由剪力產(chǎn)生的彎曲變形引起的，并不是由剪力產(chǎn)生的剪切變形引起的。在上述描述下，文獻[14]推導(dǎo)得到的法向接觸剛度系數(shù)和轉(zhuǎn)動接觸剛度系數(shù)的計算公式分別為Kθ=EI/L和Kn=EA/L，該部分的假定和推導(dǎo)是合理的，后續(xù)本文也會采用。然而，文獻[14]推導(dǎo)得到的切向彈簧剛度系數(shù)為Kτ=GA/kL，其假定和推導(dǎo)過程不合理，因此對該系數(shù)求解的修正即是本文創(chuàng)新點之一。

圖3 已有研究中的彈簧等效模型[14]Fig. 3 The spring equivalent model in existing research[14]

本文法向彈簧和轉(zhuǎn)動彈簧的定義和推導(dǎo)過程與已有研究[14]完全一致；對于切向彈簧，其變形包括由剪力產(chǎn)生的彎曲變形引起的切向位移和由剪力產(chǎn)生的剪切變形引起的附加切向位移，如圖4 所示。

圖4 等效模型中切向彈簧的定義Fig. 4 New definition of tangential spring in the equivalent model

1.1 切向接觸剛度系數(shù)推導(dǎo)

假設(shè)一個時步 Δt內(nèi)等效梁的相關(guān)物理量為：局部坐標系下接觸點C的法向位移增量 Δun和切向位移增量 Δuτ；接觸點C的轉(zhuǎn)角增量即顆粒間相對轉(zhuǎn)角 Δθ ；接觸點C的接觸力增量 ΔFn和 ΔFτ，以及接觸力矩增量 ΔM；顆粒A和顆粒B的法向位移增量和切向位移增量分別為ui和wi，i=A,B；顆粒A和顆粒B的轉(zhuǎn)角增量為 θi。根據(jù)桿系離散單元的基本假定可得接觸點C的平動位移和轉(zhuǎn)角增量與兩顆粒間的關(guān)系為：

如圖4 所示，等效模型中切向彈簧考慮了剪切變形影響，即認為切向彈簧的變形包括由剪力產(chǎn)生的彎曲變形引起的切向位移wb和由剪力產(chǎn)生的剪切變形引起的附加切向位移ws，即等效梁上任意點的切向位移w=wb+ws。此時等效梁截面轉(zhuǎn)角、剪切變形和曲率為：

切向彈簧的應(yīng)變能包括彎曲應(yīng)變能和剪切應(yīng)變能兩個部分，其表達式為：

式中，k為截面切應(yīng)力不均勻系數(shù)即截面剪切校正因子。

對式(5)的第一項即由剪力產(chǎn)生的彎曲變形在切線方向上引起的等效梁應(yīng)變能，其的求解依據(jù)桿件結(jié)構(gòu)力學(xué)中的不考慮剪切變形時彎曲梁單元相關(guān)知識[10]，則此時等效梁上任意點的切向位移即撓度函數(shù)的插值為：函數(shù)。

將式(6)和式(7)代入式(5)可得：

再由平衡條件可得，通過剪切變形得到的剪力與彎矩求導(dǎo)得到的剪力相等，即式(11)等于式(12)，得：

加之，幾何關(guān)系：

式中，b為剪切變形影響系數(shù)，取值如式(14)。

由上述推導(dǎo)過程可見，桿系離散單元法會使用到能量等效原理(即接觸單元或等效梁的應(yīng)變能等于彈簧的變形能)，在計算等效梁的應(yīng)變能時也會涉及到位移形函數(shù)，但這里的位移形函數(shù)僅用于推導(dǎo)各方向接觸剛度系數(shù)，該系數(shù)被用于計算相鄰顆粒間接觸點的接觸內(nèi)力，即位移形函數(shù)并不會在桿系離散單元計算程序中。而有限單元法引入位移形函數(shù)，是為了利用變分原理得到弱形式的等效積分方程即求解連續(xù)體的控制方程。因此，位移形函數(shù)的使用亦是桿系離散單元法與有限元法的本質(zhì)區(qū)別之一。此外，式(20)中各方向接觸剛度系數(shù)的推導(dǎo)過程確保了桿系離散元法求解桿件結(jié)構(gòu)力學(xué)行為問題的正確性和精確性。

最后，對式(20)進行擴展，得到面向空間梁單元的桿系離散單元中彈性接觸剛度系數(shù)矩陣為：

值得說明的是：本文研究對象為單層球面網(wǎng)殼結(jié)構(gòu)，該類結(jié)構(gòu)桿件的長細比較大，通常會忽略剪切變形，即后續(xù)理論和算例中均令b=0。

2 桿系離散單元精細塑性鉸模型

本文將精細塑性鉸模型引入到三維桿系離散元中，即借鑒CRC 切線模量Et的概念近似考慮由殘余應(yīng)力引起的沿桿長的逐漸屈服效應(yīng)，而截面的逐漸屈服則由拋物線函數(shù)表示。在現(xiàn)有數(shù)值方法中，最常采用精細塑性鉸模型考慮材料非線性的是非線性有限元法[8]，具體的使用步驟為：1) 通過穩(wěn)定函數(shù)考慮單元二階效應(yīng)，由此得到考慮幾何非線性的空間梁單元彈性增量力-變形關(guān)系，可見式(22)；2) 在彈性增量力-變形關(guān)系的基礎(chǔ)上引入精細化塑性鉸模型以得到空間梁單元彈塑性增量力-變形關(guān)系，可見式(23)。

借用上述非線性有限元法的思路，文中將精細塑性鉸模型引入三維桿系離散元法以考慮材料非線性。與非線性有限元法相比，桿系離散元法無需特殊處理即可考慮幾何非線性，故直接在彈性接觸本構(gòu)模型即式(21)的基礎(chǔ)上引入精細塑性鉸模型，進而建立適用于空間桿系結(jié)構(gòu)的彈塑性接觸本構(gòu)模型。對比式(21)與式(22)知，兩種方法單元剛度矩陣的表達式相似，但有以下兩點區(qū)別：1)三維桿系離散單元即式(21)無需特別考慮軸力與彎矩的耦合效應(yīng)，故在桿系離散元中認為s1y=1.0,s1z=1.0,s2y=0,s2z=0；2) 式(22)對應(yīng)的是節(jié)點變形，而式(21)對應(yīng)的是接觸點(圖2中C點)變形，故在桿系離散元中需將式(23)中的單元兩端截面剛度退化函數(shù)修正為接觸點處截面剛度退化函數(shù)。最后將上述參數(shù)取值代入式(23)并修正剛度退化函數(shù)，得到適用于空間桿系結(jié)構(gòu)的桿系離散元彈塑性接觸本構(gòu)模型的表達式為：

3 大型網(wǎng)殼結(jié)構(gòu)算例驗證與分析

3.1 K6 型單層網(wǎng)殼結(jié)構(gòu)

K6 型單層網(wǎng)殼結(jié)構(gòu)的幾何構(gòu)造如圖5 所示?？缍群褪父叻謩e為30 m 和2 m，桿件均為φ180 mm×5 mm 的圓鋼管，彈性模型E=210 GPa，屈服應(yīng)力σy=240 MPa ，剪切模量G=850 GPa，材料為理想彈塑性。底部支座為鉸接，結(jié)構(gòu)僅受頂部節(jié)點A處的豎向集中荷載作用。文獻[20]結(jié)合了塑性鉸和塑性區(qū)法的基本概念，即認為單元彈塑性僅產(chǎn)生在桿端截面，桿件其他部位仍為彈性，同時將桿端截面劃分成若干小面積模擬截面的塑性發(fā)展。文獻[16]采用桿系離散元塑性區(qū)法對該結(jié)構(gòu)進行了彈塑性分析，每個接觸截面被劃分為16 個小面積。通用有限元軟件ANSYS 中，每根桿件被劃分為4 個beam189 單元，采用弧長法進行位移追蹤，且考慮殘余應(yīng)力的影響。

圖5 K6 型單層網(wǎng)殼結(jié)構(gòu)Fig. 5 K6 single-layer reticulated shell structure

為便于文獻對比，文中桿系離散單元將結(jié)構(gòu)離散為529 個顆粒，每個桿件有4 個接觸單元，計算時步為10-4s，阻尼系數(shù)為10。采用桿系離散單元位移控制法進行加載，位移控制速度為1 mm/s。材料非線性選取Orbison[18]提出的截面極限屈服函數(shù)作為接觸截面的屈服面方程，切線模量為CRC切線模量即式(25)，截面的退化剛度系數(shù)的取值按式(26)。

桿系離散單元得到的節(jié)點A和節(jié)點B的荷載-豎向位移曲線如圖6 和圖7 所示。由圖可見，文中桿系離散單元得到的荷載-位移曲線與已有文獻、ANSYS 計算結(jié)果在彈性階段均吻合較好；相比于已有文獻[20]，桿系離散單元得到的3 條荷載-位移曲線在塑性階段與ANSYS 計算結(jié)果更吻合；當桿件進入塑性后，塑性鉸法的計算結(jié)果明顯大于塑性區(qū)法[16]和ANSYS，與兩者的最大誤差分別為9.09%和14.78%，這是因為塑性鉸法無法考慮塑性在接觸截面上和沿接觸單元(即等效梁)長的分布和發(fā)展，使得結(jié)構(gòu)的整體剛度偏大；相比于塑性鉸法，桿系離散單元精細塑性鉸法與ANSYS 計算結(jié)果和已有文獻[16]在塑性階段吻合更好，且與文獻[16]的最大誤差僅為2.2%，表明精細塑性鉸法的計算精度較塑性法有顯著提高。

圖6 節(jié)點A 的荷載-位移曲線Fig. 6 Load-displacement curve of node A

圖7 節(jié)點B 的荷載-位移曲線Fig. 7 Load-displacement curve of node B

桿系離散單元塑性鉸法、本文精細塑性鉸法以及塑性區(qū)法各自的計算耗時如表1 所示。由表可見，桿系離散單元精細塑性鉸法和塑性區(qū)法的計算時間分別是塑性鉸法的1.21 倍和2.88 倍，可見精細塑性鉸法未顯著增加桿系離散單元的計算量；與桿系離散單元塑性區(qū)法[16]相比，在誤差可接受的情況下計算耗時被大大減小，可以推測在大型結(jié)構(gòu)彈塑性分析時，計算耗時的對比會更加明顯。因此，沿桿長塑性分布不明顯時，采用桿系離散單元精細塑性鉸法求解最為劃算。

表1 桿系離散元中三種彈塑性分析方法的計算效率對比Table 1 Comparison of the calculation efficiency of three elastoplastic analysis methods using MDEM

3.2 大型單層球面網(wǎng)殼振動臺試驗

本次振動臺試驗由同濟大學(xué)多點振動臺試驗系統(tǒng)完成，該系統(tǒng)由4 個振動臺組成，試驗?zāi)Ｐ偷兹苍O(shè)40 個支座，平均落在4 個振動臺上，其余節(jié)點處于懸空狀態(tài)，如圖8 所示。試驗?zāi)Ｐ偷目s尺比為1/3.5，結(jié)構(gòu)類型為K6 型單層球面網(wǎng)殼，試驗?zāi)Ｐ蛧栏癜凑赵偷耐負潢P(guān)系，不作任何簡化。試驗?zāi)Ｐ偷目缍葹?3.4 m，矢跨比為0.5，焊接空心球數(shù)為1261，桿件數(shù)為3660，模型節(jié)點的配重為30 kg。桿件截面共有4 種規(guī)格，且滿足滿應(yīng)力設(shè)計準則，其中，截面尺寸為φ23 mm×1 mm 、φ38 mm×2 mm 的桿件為采用經(jīng)過真空退火處理的20#鋼，以模擬原型材料Q235 的應(yīng)力-應(yīng)變關(guān)系，其余桿件采用Q235 熱軋鋼管，每種規(guī)格桿件的各項力學(xué)指標詳見文獻[21 - 22]。

圖8 振動臺布置與試驗?zāi)Ｐ虵ig. 8 Shaking table layout and test model

試驗采用EI-Centro 地震波(1940)，根據(jù)模型相似比將該波的持時壓縮為28 s，時間間隔為0.010 69 s。每個振動臺實際輸入的地震波加速度時程曲線和速度時程曲線、以及加載工況詳見文獻[21 - 22]。

試驗?zāi)Ｐ捅浑x散為14 581 個顆粒和16 980 個接觸單元，得到的桿系離散模型如圖9 所示。采用位移法施加多點激勵的地震作用，具體加載的步驟和公式見文獻[22]。結(jié)構(gòu)的阻尼比和基頻根據(jù)白噪聲掃頻結(jié)果計算得到，數(shù)值如圖10 所示。材料本構(gòu)為理想彈塑性模型，仍選取Orbison[18]提出的截面極限屈服函數(shù)作為接觸截面的屈服面方程。

圖9 試驗?zāi)Ｐ偷臈U系離散單元數(shù)值模型Fig. 9 Discrete element model of test model

圖10 試驗?zāi)Ｐ偷淖枘岜群突lFig. 10 Damping ratio and fundamental frequency of test model

圖11 給出了試驗?zāi)Ｐ偷腜GA-Δmax曲線對比(荷載-位移選取點見圖9)，由圖可見，數(shù)值仿真結(jié)果和試驗結(jié)果的變化趨勢相同，即隨著峰值加速度的提高，最大位移顯著增加，相應(yīng)倒塌的前一級加載工況，模型節(jié)點最大位移達到70.71 mm(試驗為78.76 mm)，為跨度的1 /330，倒塌征兆顯著，PGA-Δmax曲線斜率逐漸減小，未有突變現(xiàn)象，結(jié)構(gòu)表現(xiàn)出良好的延性，塑性發(fā)展充分，說明結(jié)構(gòu)的破壞屬于承載力破壞。此外，從圖11 的對比結(jié)果亦可得出與算例3.1 一致的結(jié)論，即當結(jié)構(gòu)進入動力彈塑性行為后，桿系離散元精細塑性鉸法的計算精度明顯高于塑性鉸法，進一步驗證了前者的可靠性。

圖11 PGA-Δmax 曲線對比Fig. 11 Comparison of PGA-Δmax curves

將試驗?zāi)Ｐ驮诓煌r下位移時程曲線與試驗結(jié)果進行對比如圖12～圖13。由圖可見，桿系離散單元精細塑性鉸法計算得到的位移時程反應(yīng)曲線與試驗結(jié)果幾乎重合，僅在峰值上有些許誤差，可能原因是：1)試驗中各工況下使用模型的初始狀態(tài)為前一級工況下變形后的試驗?zāi)Ｐ?，而?shù)值仿真中各工況下采用的計算結(jié)構(gòu)均為未變形的原始模型；2)數(shù)值仿真中未考慮地震荷載下鋼材應(yīng)變率效應(yīng)的影響，使得部分桿件的屈服強度低于實際值，使得數(shù)值計算位移大于試驗值，這一現(xiàn)象在變形大的工況較為明顯，如圖12(b)、圖12(c)所示。

圖12 不同工況X 向最大位移測點的時程曲線位移對比Fig. 12 Comparison of displacement time history curves of the maximum displacement of measuring point in X-direction in different cases

圖13 不同工況Y 向最大位移測點的時程曲線位移對比Fig. 13 Comparison of displacement time history curves of the maximum displacement of measuring point in Y-direction in different cases

圖14 給出了400 gal 工況下測點的位移時程曲線。由圖可知：數(shù)值仿真得到的模型極限荷載為400 gal，與試驗結(jié)果一致，但二者完全倒塌時間點略有差別，試驗中完全倒塌時刻為6.62 s，數(shù)值模擬中為5.69 s，誤差較小。

圖14 400 gal 工況下測點的位移時程曲線Fig. 14 Displacement time history curve of a measuring point in the case of (PGA=400 gal)

圖15 為400 gal 極限荷載下試驗?zāi)Ｐ推茐倪^程的數(shù)值仿真結(jié)果，該仿真結(jié)果與試驗結(jié)果[21-22]相一致。由圖可知：隨著輸入峰值加速度的提高，各臺面上方第2 圈、第3 圈斜桿首先發(fā)生彎曲；之后，結(jié)構(gòu)整體變形程度逐漸加大、破壞區(qū)域逐漸擴展，同時與地震波輸入方向垂直的臺面之間第3 圈～第5 圈桿件變形明顯，但上部仍處于彈性狀態(tài)，無明顯變形；最終，模型支座上方及與地震波輸入方向垂直的臺面之間第1 圈～第6 圈桿件均變形顯著，整個結(jié)構(gòu)向下坍塌。可見，試驗?zāi)Ｐ偷臈U件屈曲過程是由下而上、由各臺面上方破壞區(qū)域向兩側(cè)擴展；對于正常設(shè)計的網(wǎng)殼結(jié)構(gòu)，支座上方桿件和與地震波輸入方向垂直的主肋區(qū)下部桿件及其附近區(qū)域為結(jié)構(gòu)關(guān)鍵部位。

圖15 試驗?zāi)Ｐ偷钠茐倪^程數(shù)值仿真結(jié)果Fig. 15 Numerical simulation results of the failure process of the test model

4 結(jié)論

本文對桿系離散單元法的基本計算理論進行了修正，并將精細塑性鉸法引入其中，最后通過兩個大型網(wǎng)殼結(jié)構(gòu)的彈塑性行為驗證了所提方法的精確性，得到如下結(jié)論：

(1)重新定義了桿系離散單元中接觸本構(gòu)模型的切向彈簧，并根據(jù)能量等效原理推導(dǎo)了不考慮剪切變形和考慮剪切變形時切向接觸剛度系數(shù)的統(tǒng)一計算公式，實現(xiàn)了桿系離散單元基本理論系統(tǒng)化，為結(jié)構(gòu)力學(xué)行為模擬提供更為嚴謹?shù)睦碚撝巍?/p>

(2)提出了桿系離散單元精細塑性鉸法，并推導(dǎo)了精細塑性鉸法的彈塑性接觸本構(gòu)模型。其通過切線模量和截面剛度退化系數(shù)近似考慮殘余應(yīng)力對接觸單元剛度的削弱，前者考慮了沿接觸單元長的逐漸屈服效應(yīng)；后者表示接觸截面的逐漸屈服。靜力彈塑性分析算例表明，桿系離散單元精細塑性鉸法可近似考慮構(gòu)件的塑性發(fā)展，其計算精度明顯高于塑性鉸法，且不會顯著增加桿系離散單元的計算量；當截面分布塑性不明顯時，相比于塑性區(qū)法，采用桿系離散單元精細塑性鉸法“性價比”更高。

(3)將修正后的桿系離散單元法應(yīng)用于多點激勵單層球面網(wǎng)殼彈塑性行為模擬，并與試驗結(jié)果對比，分析結(jié)果表明，二者極限荷載、倒塌模式一致，且數(shù)值仿真得到的位移時程曲線與試驗值幾乎重合，僅在幅值上有微小差異，進一步驗證了本文方法處理大跨空間結(jié)構(gòu)動力彈塑性問題的精確性。