柔性框架結(jié)構(gòu)動(dòng)力非線性分析的剛體準(zhǔn)則法

陳朝暉，陶宇宸，何 敏

(1. 重慶大學(xué)土木工程學(xué)院，重慶 400045；2. 山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))，重慶 400045；3. 浙江大學(xué)建筑工程學(xué)院，浙江，杭州 310058)

高層、高聳、大跨結(jié)構(gòu)以及空間網(wǎng)殼結(jié)構(gòu)等柔性結(jié)構(gòu)，在地震、風(fēng)、海浪等動(dòng)力荷載作用下的大位移、大轉(zhuǎn)動(dòng)等幾何非線性特征顯著，且?guī)缀畏蔷€性與動(dòng)力效應(yīng)相互耦合，給精確而高效的數(shù)值分析造成了困難。

結(jié)構(gòu)運(yùn)動(dòng)方程的求解通常采用直接積分法，包括顯式積分法和隱式積分法。顯式算法較為簡單，但對于強(qiáng)非線性問題，需采用較小的時(shí)間步長來保證計(jì)算精度，計(jì)算效率低。隱式算法因在對未知時(shí)間步的求解中，包含了與該時(shí)間步相關(guān)的一個(gè)或多個(gè)未知量，因而需要迭代求解。對于強(qiáng)非線性問題，長時(shí)間的響應(yīng)分析會(huì)大到不切實(shí)際[1]。非線性動(dòng)力分析的關(guān)鍵之一是如何將單元非線性位形描述與運(yùn)動(dòng)方程的求解結(jié)合。常用方法包括基于TL 列式、CR 列式以及二者混合的方法。Bathe 等[2]最早引入TL 列式的靜力非線性分析方法，建立了動(dòng)力大變形問題的有限元分析方法。在此基礎(chǔ)上，Remseth[3]采用TL 列式分析了結(jié)構(gòu)幾何非線性動(dòng)力響應(yīng)。但TL 列式始終以初始構(gòu)型為參考，當(dāng)柔性結(jié)構(gòu)具有較大的振動(dòng)位移時(shí)，計(jì)算結(jié)果與實(shí)際情況偏差較大。為此，有研究者[4]提出了CR 列式的非線性動(dòng)力分析方法，將單元變形分解為剛體位移和自然變形，通過扣除初始狀態(tài)到當(dāng)前位形的剛體位移來得到單元結(jié)點(diǎn)實(shí)際位移。但這類CR 列式法在工程應(yīng)用中遭遇很大局限，原因在于剛體位移和自然變形的分解導(dǎo)致單元?jiǎng)恿?xiàng)的表達(dá)式推導(dǎo)非常復(fù)雜，難以為工程接受。有研究者提出了CR 列式與TL 列式結(jié)合的方法，如Le 等[5-6]利用CR 列式推導(dǎo)了單元的剛度矩陣，同時(shí)使用TL 列式推導(dǎo)了慣性力向量，建立了平面梁和空間梁的動(dòng)力分析模型。

幾何非線性問題的難點(diǎn)在于單元的大變形和大轉(zhuǎn)動(dòng)會(huì)造成附加內(nèi)力，若不能合理描述單元的變形及其產(chǎn)生的結(jié)點(diǎn)力增量，其誤差經(jīng)累計(jì)后將使計(jì)算結(jié)果嚴(yán)重偏離實(shí)際。近年來，基于CR 列式的非線性分析方法發(fā)展迅速，并用于彈塑性分析中[7-8]。通常，基于CR 列式的單元著眼于對單元變形的描述，導(dǎo)致變形描述準(zhǔn)確的單元過于復(fù)雜，而變形近似的單元精度又差強(qiáng)人意。事實(shí)上，對于初始平衡的單元，若僅發(fā)生剛體轉(zhuǎn)動(dòng)，其平衡的結(jié)點(diǎn)內(nèi)力必將隨單元發(fā)生剛體轉(zhuǎn)動(dòng)，而大小不變，從而使單元在當(dāng)前狀態(tài)下繼續(xù)維持平衡。此即大位移大轉(zhuǎn)動(dòng)分析的“剛體準(zhǔn)則”，由Yang 等[9]于1987 年率先提出。基于剛體準(zhǔn)則，筆者先后建立了一系列線彈性桁架單元[10]、平面與空間梁單元[10-11]以及板、殼及膜單元[12-13]。并結(jié)合塑性鉸理論，將剛體準(zhǔn)則推廣至柔性框架結(jié)構(gòu)的彈塑性非線性靜力分析[14-15]。上述滿足剛體準(zhǔn)則的各類單元與其他分析方法與商業(yè)軟件相比，精度與效率優(yōu)勢顯著。

鑒于上述剛體準(zhǔn)則及其相應(yīng)單元在靜力幾何非線性分析中的優(yōu)勢，本文提出了一種高效且高精度的柔性空間桿系結(jié)構(gòu)動(dòng)力非線性分析方法。該方法采用滿足剛體準(zhǔn)則的空間梁單元，采用HHT-α 隱式積分法將運(yùn)動(dòng)方程轉(zhuǎn)化為等效動(dòng)力增量方程，進(jìn)而將剛體準(zhǔn)則融入求解等效增量方程的Newton-Raphson 方法，從而建立了柔性框架結(jié)構(gòu)的動(dòng)力時(shí)程分析方法。典型柔性框架結(jié)構(gòu)動(dòng)力分析及其與Le 等[5-6]、Cho 等[16]高精度方法，以及ABAQUS 商業(yè)軟件等的對比表明，本文方法對大位移下的結(jié)點(diǎn)力增量計(jì)算簡潔，單元數(shù)和迭代步少，精度高，適于工程應(yīng)用。

1 等效動(dòng)力增量方程

結(jié)構(gòu)系統(tǒng)運(yùn)動(dòng)方程的一般形式可記作：

式中：M、C和K分別為結(jié)構(gòu)整體質(zhì)量陣、阻尼陣和剛度陣；Pt為動(dòng)力荷載向量；Ut、U˙t和U¨t為位移向量、速度向量和加速度向量。式中右上標(biāo)t為時(shí)間，下文同。

直接積分法是求解動(dòng)力問題的常用方法，其基本思想基于差分法，即將計(jì)算總時(shí)長劃分為若干時(shí)間步t1,t2, ···,tn，在每一時(shí)間步Δt=ti-ti-1內(nèi)，人為假設(shè)位移u、速度u˙ 和加速度u¨的關(guān)系，而運(yùn)動(dòng)方程僅在各時(shí)間步的兩端滿足。由Hilber 和Hughes 提出的HHT-α 法[17]通過參數(shù) α來引入數(shù)值阻尼，從而防止計(jì)算發(fā)散，選取相對較大步長時(shí)仍能保持?jǐn)?shù)值穩(wěn)定。設(shè)每一時(shí)間步內(nèi)，位移、速度與加速度的近似關(guān)系為：

式(9)形同增量形式的靜力平衡方程，因而也稱為擬靜力增量方程。籍此，就將求解運(yùn)動(dòng)方程的積分問題轉(zhuǎn)化為在每一個(gè)時(shí)間步內(nèi)求解非線性靜力增量方程的問題。

考慮到阻尼的影響，還需引入阻尼模型。在此，采用結(jié)構(gòu)工程中較通行的Rayleigh 阻尼模型，即：

式中，a0和a1為比例常數(shù)。

由于直接積分法基于差分法，運(yùn)動(dòng)方程只在每個(gè)時(shí)間步的頭、尾兩端嚴(yán)格滿足，而在時(shí)間步內(nèi)存在人為近似，從而導(dǎo)致求解過程中的非線性，同時(shí)引入阻尼會(huì)帶來動(dòng)力響應(yīng)的非線性特性。換言之，即使對于M和K為常數(shù)陣的動(dòng)力問題，在采用直接積分法并考慮阻尼的情況下，系統(tǒng)的動(dòng)力響應(yīng)也表現(xiàn)出非線性特性。對于具有大位移和大轉(zhuǎn)動(dòng)的柔性結(jié)構(gòu)，其動(dòng)力響應(yīng)的非線性還體現(xiàn)在其他兩方面。其一，是結(jié)構(gòu)整體剛度矩陣K依賴于結(jié)構(gòu)的瞬時(shí)位形，K不再是常數(shù)而隨位形變化，剛度矩陣的非線性同時(shí)造成阻尼矩陣C的非線性。若采用一致質(zhì)量模型，則結(jié)構(gòu)質(zhì)量矩陣M也是非線性的。本文為簡化起見，采用集中質(zhì)量模型。其二，與靜力非線性分析類似，結(jié)構(gòu)動(dòng)力響應(yīng)中的大位移和大轉(zhuǎn)動(dòng)還會(huì)引起附加結(jié)點(diǎn)力，該附加結(jié)點(diǎn)力若不能合理計(jì)算，將使系統(tǒng)等效動(dòng)力增量方程在每一個(gè)時(shí)間步內(nèi)無法平衡，而該不平衡力將無法通過減小時(shí)間步長來消除。顯然，大位移大轉(zhuǎn)動(dòng)情形下的結(jié)構(gòu)非線性動(dòng)力響應(yīng)問題求解的關(guān)鍵在于結(jié)構(gòu)切線剛度陣的確定以及結(jié)點(diǎn)力增量的計(jì)算。由式(9)增量形式的非線性動(dòng)力平衡方程可以看出，在UL 列式下，該非線性動(dòng)力問題求解過程中的結(jié)點(diǎn)力增量計(jì)算與靜力非線性問題本質(zhì)是相同的，因此，接下來即討論如何將靜力非線性分析的剛體準(zhǔn)則植入直接積分法中，從而建立動(dòng)力非線性問題的迭代求解方法。

2 動(dòng)力非線性分析的剛體準(zhǔn)則法

2.1 單元結(jié)點(diǎn)力增量的計(jì)算

如圖1 所示的壓桿屈曲問題，可將單元的大變形和大轉(zhuǎn)動(dòng)視為兩個(gè)過程的組合[10]：單元先發(fā)生由初始平衡狀態(tài)C1至當(dāng)前狀態(tài)C2的剛體轉(zhuǎn)動(dòng)ur，而后在C2狀態(tài)下產(chǎn)生彈性變形，稱其為“自然變形”un，則單元變形可寫作：

圖1 懸臂壓桿屈曲變形[11]Fig. 1 Deformation of a buckling cantilever

在剛體轉(zhuǎn)動(dòng)階段，C1狀態(tài)的單元結(jié)點(diǎn)力僅隨單元發(fā)生剛體位移，大小不變，在C2狀態(tài)下仍然平衡，如圖2 所示空間梁單元；單元的結(jié)點(diǎn)力增量由C2狀態(tài)下的自然變形un產(chǎn)生。

圖2 經(jīng)歷剛體轉(zhuǎn)動(dòng)的空間梁單元[11]Fig. 2 Three-dimensional beam element experiencing rigid rotation

不失合理性，可以認(rèn)為對于大多數(shù)工程大位移大轉(zhuǎn)動(dòng)問題，剛體位移占單元位移的主要部分，相較之下當(dāng)前狀態(tài)的自然變形是小量。對于線彈性問題，基于剛體準(zhǔn)則，單元結(jié)點(diǎn)力在剛體位移上不做功，即keur=0(ke為單元彈性剛度矩陣)，單元的結(jié)點(diǎn)力增量 Δf僅由單元的線彈性變形引起，為：

式(16)表示，將上一狀態(tài)平衡的單元結(jié)點(diǎn)力轉(zhuǎn)動(dòng)至當(dāng)前位置1f，再疊加彈性變形引起的結(jié)點(diǎn)力增量 Δf。

進(jìn)一步地，可由虛功原理建立單元增量平衡方程[7]：

式中，kg為單元幾何剛度矩陣，具體形式參見文獻(xiàn)[11]。

對于伴隨大位移大轉(zhuǎn)動(dòng)的結(jié)構(gòu)非線性動(dòng)力問題，可以利用上述剛體準(zhǔn)則來處理每一時(shí)間步內(nèi)的結(jié)點(diǎn)力增量，即在每一時(shí)間步內(nèi)，可認(rèn)為單元的剛體轉(zhuǎn)動(dòng)在其位移增量中占比較大，而自然變形相對較小，材料線性條件下則為線彈性小變形。因此，可將上一時(shí)間步末滿足平衡條件的單元結(jié)點(diǎn)力保持其大小不變，而隨單元旋轉(zhuǎn)至當(dāng)前狀態(tài)，再疊加由單元彈性變形引起的結(jié)點(diǎn)力增量。當(dāng)前時(shí)間步下單元的位置則由HHT-α 法確定。

2.2 動(dòng)力非線性剛體準(zhǔn)則法基本步驟

圖3 動(dòng)力非線性分析的擬靜力增量-迭代法示意圖Fig. 3 Scheme of quasi static incremental-iteration method for nonlinear dynamic analysis

圖4 給出了上述植入剛體準(zhǔn)則的空間柔性桿系結(jié)構(gòu)動(dòng)力非線性分析流程。

圖4 動(dòng)力非線性分析計(jì)算流程圖Fig. 4 Flow graph of nonlinear dynamic analysis method

3 算例與分析

本節(jié)通過2 個(gè)柔性框架結(jié)構(gòu)算例驗(yàn)證本文方法對于線彈性幾何非線性動(dòng)力問題的有效性，并與已有文獻(xiàn)和ABAQUS 結(jié)果對比，驗(yàn)證本文方法在動(dòng)力非線性問題上的精度與效率。各算例分析中HHT-α 法中的參數(shù) Δt取為-0.01，在計(jì)算過程引入微小數(shù)值阻尼過濾高頻響應(yīng)，但是數(shù)值阻尼會(huì)引起能量耗散，在時(shí)間步較大時(shí)會(huì)影響結(jié)果精度，所以在確定各算例時(shí)間步長 Δt時(shí)都采用較小步長 0.5Δt進(jìn)行驗(yàn)證，若兩次得到的結(jié)果相同即可確定步長為 Δt。

3.1 空間懸臂梁

圖5 所示為自由端受簡諧荷載激勵(lì)的懸臂梁，梁的跨度L=10 m，荷載幅值為10 MN，頻率為50 rad/s，梁截面尺寸為0.25 m ×0.5 m，彈性模量E=210 GPa，材料密度 ρ=7850 kg/m3。該懸臂梁是檢驗(yàn)幾何非線性方法的經(jīng)典算例，Cho 等[16]和Le 等[5]采用CR 列式進(jìn)行了位移時(shí)程分析。在此，本文基于前述剛體準(zhǔn)則非線性動(dòng)力分析方法，采用文獻(xiàn)[18]所建滿足剛體準(zhǔn)則的歐拉梁單元，分析該懸臂梁在簡諧激勵(lì)下的動(dòng)力響應(yīng)。時(shí)間步長 Δt取為10-5s，劃分為3 個(gè)單元。同時(shí)與Cho 等[16]以及ABAQUS 軟件計(jì)算結(jié)果進(jìn)行對比，其中ABAQUS 采用Beam21 平面梁單元，劃分10 個(gè)單元，同樣采用HHT-α 法計(jì)算。

圖5 懸臂梁幾何尺寸Fig. 5 Geometrical data of cantilever beam

圖6 和圖7 為懸臂梁自由端豎向與水平位移響應(yīng)時(shí)程，可以看出，梁自由端振幅較大，幾何非線性特征明顯。由圖5 可以看出，采用ABAQUS分析時(shí)考慮幾何非線性與否，結(jié)果差異顯著。本文方法與Cho 等[16]、ABAQUS 非線性分析結(jié)果在數(shù)值及變化趨勢上完全一致。ABAQUS 用了10 單元，而本文只需要3 個(gè)單元；Cho 等[16]采用CR列式推導(dǎo)的三角形單元，將結(jié)構(gòu)劃分為24 個(gè)平面三角形單元分析，單元復(fù)雜，計(jì)算成本高。因此，本文方法對于大轉(zhuǎn)動(dòng)大位移的非線性動(dòng)力問題，在計(jì)算精度和效率兩方面均具有顯著優(yōu)勢。

圖6 懸臂梁豎向位移時(shí)程曲線Fig. 6 Vertical displacement history of cantilever beam

圖7 懸臂梁水平位移時(shí)程曲線Fig. 7 Horizontal displacement history of cantilever beam

3.2 Lee 框架

如圖8 所示，由2 根正交矩形截面直桿組成的兩鉸框架，水平桿在距兩桿交點(diǎn)L/5 處受豎直向下的集中荷載作用，此框架結(jié)構(gòu)被稱為“Lee 框架”。由于該框架表現(xiàn)出的復(fù)雜后屈曲行為，也被視作檢驗(yàn)非線性方法合理性的典型問題。Le 等[5]采用CR 列式對此進(jìn)行了突加荷載下的動(dòng)力響應(yīng)分析，每根桿件劃分了5 個(gè)單元。Lee 框架每根桿長L=12 m，梁橫截面尺寸為a=0.2 m、e=0.3 m，彈性模量E=210 GPa，密度為 ρ=7850 kg/m3。水平桿A 點(diǎn)作用一豎直向下的突加荷載P，大小為4.1 MN。在此，采用本文所提出的基于剛體準(zhǔn)則的動(dòng)力非線性分析方法，每根桿件劃分為10 個(gè)單元，時(shí)間步長 Δt取為5×10-5s，分析荷載作用點(diǎn)A的豎向和水平位移時(shí)程。作為對比的ABAQUS分析采用Beam21 梁單元，每根桿件等分為10 個(gè)單元。

圖8 Lee 框架Fig. 8 Geometrical data of Lee's frame

圖9 和圖10 分別為Lee 框架A點(diǎn)的豎向與水平位移時(shí)程，在突加荷載后的前1.5 s 內(nèi)，三種方法結(jié)果基本吻合，而在荷載作用后期，本文方法與ABAQUS 的結(jié)果仍然吻合較好，TN-Le 的結(jié)果與ABAQUS 結(jié)果偏差較大。

圖9 突加荷載作用下Lee 框架A 點(diǎn)豎向位移Fig. 9 Vertical displacement history of A point of Lee's frame under sudden load

圖10 突加荷載作用下Lee 框架A 點(diǎn)水平位移Fig. 10 Horizontal displacement history of A point of Lee's frame under sudden load

進(jìn)一步地，將荷載更換為A點(diǎn)施加豎直向下的簡諧荷載P=4.1×106sin(50t)，計(jì)算得到A點(diǎn)的豎向與水平位移時(shí)程曲線如圖11～圖12 所示。可以看出本文方法與ABAQUS 每根桿劃分10 個(gè)單元的結(jié)果一致，表明本文方法對不同荷載工況下的動(dòng)力響應(yīng)分析問題均能適用。

圖11 簡諧荷載作用下Lee 框架A 點(diǎn)豎向位移Fig. 11 Vertical displacement history of A point of Lee's frame under harmonic load

圖12 簡諧荷載作用下Lee 框架A 點(diǎn)水平位移Fig. 12 Horizontal displacement history of A point of Lee's frame under harmonic load

4 結(jié)論

本文基于靜力幾何非線性分析的剛體準(zhǔn)則，結(jié)合HHT- α 隱式積分法，推導(dǎo)了動(dòng)力分析的有限元增量求解格式，建立了簡潔高效的柔性框架結(jié)構(gòu)動(dòng)力響應(yīng)分析方法，主要結(jié)論如下：

(1)為了考慮空間柔性桿系結(jié)構(gòu)大幅振動(dòng)的幾何非線性效應(yīng)，本文在動(dòng)力時(shí)程計(jì)算的每一個(gè)時(shí)間步中，考慮結(jié)構(gòu)的大位移大轉(zhuǎn)動(dòng)，在每一步的迭代中采用剛體準(zhǔn)則來處理?xiàng)U件的剛體轉(zhuǎn)動(dòng)。該方法可以有效分析柔性空間框架結(jié)構(gòu)的幾何非線性動(dòng)力問題，合理描述結(jié)構(gòu)在簡諧荷載、突加荷載等動(dòng)力荷載作用下的結(jié)構(gòu)響應(yīng)。

(2)在計(jì)算單元結(jié)點(diǎn)力時(shí)，將單元變形過程看作由初始平衡位置發(fā)生的剛體位移以及在運(yùn)動(dòng)后位置上的自然變形，初始單元平衡力隨著剛體位移的過程發(fā)生隨動(dòng)，而大小不變，在此基礎(chǔ)上疊加自然變形引起的結(jié)點(diǎn)力增量，計(jì)算精度與效率高。

(3)通過分析具有大變形大轉(zhuǎn)角等典型幾何非線性特征的桿系結(jié)構(gòu)，驗(yàn)證了本文推導(dǎo)的剛體準(zhǔn)則單元對于動(dòng)力大變形問題的適用性，具有極高的分析精度，克服了以往方法對于轉(zhuǎn)角限制的假設(shè)。柔性空間桿系結(jié)構(gòu)的算例分析，體現(xiàn)了本文方法在計(jì)算效率上相對于現(xiàn)有研究方法以及有限元軟件的優(yōu)勢，單元?jiǎng)澐稚伲啾扔趥鹘y(tǒng)的TL 列式、CR 列式計(jì)算效率大大提升。對于不同的工程結(jié)構(gòu)，只需修改結(jié)構(gòu)模型以及力和位移的邊界條件，易于程序編制，避免類似于CR 列式對于結(jié)點(diǎn)轉(zhuǎn)動(dòng)向量的儲(chǔ)存，迭代收斂快。