柔性框架結(jié)構動力非線性分析的剛體準則法

陳朝暉，陶宇宸，何 敏

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

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

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

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

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

1 等效動力增量方程

結(jié)構系統(tǒng)運動方程的一般形式可記作：

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

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

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

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

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

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

2 動力非線性分析的剛體準則法

2.1 單元結(jié)點力增量的計算

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

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

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

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

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

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

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

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

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

2.2 動力非線性剛體準則法基本步驟

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

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

圖4 動力非線性分析計算流程圖Fig. 4 Flow graph of nonlinear dynamic analysis method

3 算例與分析

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

3.1 空間懸臂梁

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

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

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

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

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

3.2 Lee 框架

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

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

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

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

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

進一步地，將荷載更換為A點施加豎直向下的簡諧荷載P=4.1×106sin(50t)，計算得到A點的豎向與水平位移時程曲線如圖11～圖12 所示?？梢钥闯霰疚姆椒ㄅcABAQUS 每根桿劃分10 個單元的結(jié)果一致，表明本文方法對不同荷載工況下的動力響應分析問題均能適用。

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

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

4 結(jié)論

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

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

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

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