謝金哲， 吳 斌,2， 楊浩文

(1. 哈爾濱工業(yè)大學(xué) 土木工程學(xué)院, 哈爾濱 150090； 2. 武漢理工大學(xué) 土木工程與建筑學(xué)院, 武漢 430070)

逐步積分算法是求解結(jié)構(gòu)動(dòng)力響應(yīng)的重要手段之一，它同時(shí)也被運(yùn)用到混合試驗(yàn)這類試驗(yàn)與數(shù)值計(jì)算相結(jié)合的試驗(yàn)中[1-2]。而無條件穩(wěn)定的逐步積分算法能夠保證混合試驗(yàn)這類極度依賴算法穩(wěn)定性的試驗(yàn)順利進(jìn)行，因此在此類依賴算法穩(wěn)定性的計(jì)算中具有巨大優(yōu)勢。

最常用的逐步積分算法為中心差分法和Newmark-β[3]方法，其穩(wěn)定性在線性體系下已經(jīng)得到了充分的研究。但是從理論上講，諸如譜半徑方法[4]和極點(diǎn)法[5]等線性系統(tǒng)數(shù)值穩(wěn)定性判定方法在非線性系統(tǒng)下并不適用。Xie等[6]通過算例指出線性體系下無條件穩(wěn)定的算法，比如平均加速度方法，在非線性領(lǐng)域并不一定具有無條件穩(wěn)定性。Belytschko等[7]指出離散系統(tǒng)的能量有界是算法保持穩(wěn)定性的充分條件，此判定準(zhǔn)則對非線性系統(tǒng)也適用。潘天林等[8]采用能量守恒準(zhǔn)則討論了隱式中點(diǎn)法對非線性指數(shù)阻尼的結(jié)構(gòu)動(dòng)力方程為數(shù)值穩(wěn)定。但是Belytschko等把能量有界作為現(xiàn)有方法的方程組迭代收斂收斂條件，可能會(huì)導(dǎo)致方程組無解。此后的學(xué)者將能量守恒或能量一致作為約束條件，同時(shí)引入新的參數(shù)對現(xiàn)有算法進(jìn)行修正，為新算法的可解性提供了可能性。

最早的能量積分方法由Hughes等[9]提出，在平均加速度方法的基礎(chǔ)之上，將能量守恒視為一種約束條件，將平衡方程視為某種泛函的條件極值，通過找到滿足變分關(guān)系的泛函表達(dá)式，將時(shí)間離散問題轉(zhuǎn)化為泛函的條件極值問題，使用拉格朗日乘子法對平衡方程進(jìn)行修正。吳斌等[10]討論了保留平衡方程，修改速度-加速度關(guān)系的Hughes方法，但是這樣做會(huì)導(dǎo)致局部精度可能退化到0階。Crisfield等[11]在94年將共旋方法運(yùn)用到桁架體系并提出了考慮幾何非線性的線彈性桁架的能量守恒方法，該方法只能運(yùn)用在線彈性體系。Simo等[12-13]通過對隱式中點(diǎn)法進(jìn)行修正實(shí)現(xiàn)了三維實(shí)體單元和幾何精確梁單元的具有二階局部精度的能量守恒積分方法,張閏等[14]將求積元運(yùn)用到Simo提出的幾何精確梁能量方法中。李妍等[15]對比了Simo方法和Hughes方法，理論分析與數(shù)值算例結(jié)果表明Simo能量法優(yōu)于Hughes能量法和平均加速度法。幾何精確梁在有限元計(jì)算的時(shí)候通過對轉(zhuǎn)角和軸線的獨(dú)立插值，在運(yùn)用能量方法時(shí)可以得到簡潔的表達(dá)式，但是對于轉(zhuǎn)角和位移具有微分關(guān)系的歐拉梁而言，直接套用幾何精確梁的方法是不合適的。Galvanetto等[16]通過對隱式中點(diǎn)法進(jìn)行修正得到了共旋歐拉梁的能量守恒方法，為了得到修正參數(shù)的顯式表達(dá)式，他們對剛體轉(zhuǎn)角的變化量進(jìn)行了近似處理，當(dāng)軸向變形較大時(shí)候該算法無法保證能量守恒。潘天林[17]在文獻(xiàn)[13]的基礎(chǔ)之上提出了平面共旋歐拉梁的非線性修正方法，通過對幾何非線性和材料非線性項(xiàng)分別進(jìn)行修正解決了文獻(xiàn)[13]的問題，給出了具有二階局部精度的能量方法。Chhang等[18]提出了具有二階局部精度的線彈性平面共旋歐拉梁的能量-動(dòng)量守恒方法，為了使算法滿足能量守恒，他用平均應(yīng)變近似中點(diǎn)構(gòu)型下的應(yīng)變，本質(zhì)上是對平衡方程的修正。目前尚未見到大變形下歐拉梁單元能量守恒數(shù)值算法方面的研究。

本文根據(jù)文獻(xiàn)[19]提出的對水平和豎向位移場獨(dú)立插值方法，考慮大轉(zhuǎn)角的影響并用工程應(yīng)變描述軸向應(yīng)變，實(shí)現(xiàn)對大變形歐拉梁的描述。在時(shí)間離散上采用文獻(xiàn)[17]提出的非線性修正方法。但是本文提出的空間離散方法中材料非線性項(xiàng)和幾何非線性項(xiàng)無法分離，因此無法實(shí)現(xiàn)文獻(xiàn)[17]中單元層面上的對幾何非線性項(xiàng)和材料非線性項(xiàng)分別修正的多參數(shù)修正，本文將對整體結(jié)構(gòu)采用單參數(shù)修正確保離散系統(tǒng)的能量守恒。

本文研究內(nèi)容如下：推導(dǎo)用于描述大變形歐拉梁的有限元模型，根據(jù)文獻(xiàn)[17]給出的非線性修正方法推導(dǎo)出用于本有限元模型的逐步積分方法，最后通過三個(gè)算例，驗(yàn)證該算法的有效性。

1 能量方法推導(dǎo)

1.1 空間離散

假定梁的初始構(gòu)型和當(dāng)前構(gòu)型，如圖1所示。

引入一個(gè)無量綱參數(shù)s∈[0,1]，它表示每個(gè)點(diǎn)對應(yīng)在梁上的位置，s(0)代表梁左端點(diǎn)，s(1)代表梁右端點(diǎn)，令初始構(gòu)型梁上點(diǎn)的坐標(biāo)與參數(shù)s滿足函數(shù)關(guān)系：

圖1 初始構(gòu)型與當(dāng)前構(gòu)型圖

(1)

式中：下標(biāo)o表示初始構(gòu)型，n表示當(dāng)前構(gòu)型，l為單元初始長度，梁單元的位移場滿足式(2)所示函數(shù)關(guān)系：

(2)

式中：參數(shù)a為待定參數(shù)，通過代入到邊界條件得到變形后的構(gòu)型如下：

(3)

由于在運(yùn)用虛位移方程的時(shí)候，本文采用基于初始構(gòu)型的虛位移方程，因此曲率的表達(dá)式為：

(4)

軸向應(yīng)變的表達(dá)式為：

(5)

通過對初始構(gòu)型使用虛位移原理得到結(jié)構(gòu)反力表達(dá)式：

(6)

對于線彈性本構(gòu)，其結(jié)構(gòu)反力表達(dá)式為：

(7)

式中：E為單元彈性模量，I為截面慣性矩，A為單元截面面積。由此得到梁單元的等效節(jié)點(diǎn)反力與節(jié)點(diǎn)位移之間的函數(shù)關(guān)系。同時(shí)為了運(yùn)用能量方法，需要給出計(jì)算單元的變形能的計(jì)算表達(dá)式。對于矩形截面，變形能EP與節(jié)點(diǎn)位移的函數(shù)關(guān)系為：

(8)

式中：b為截面寬度；h為截面高度；E為單元彈性勢能。

上述推導(dǎo)是在局部坐標(biāo)進(jìn)行的，對于多單元的整體坐標(biāo)，可以通過如下轉(zhuǎn)換得到整體坐標(biāo)下的位移對應(yīng)的結(jié)構(gòu)反力。整體坐標(biāo)下單元節(jié)點(diǎn)位移uG與局部坐標(biāo)單元節(jié)點(diǎn)位移uL存在轉(zhuǎn)換關(guān)系：

uL=TuG

(9)

式中：T為坐標(biāo)轉(zhuǎn)換矩陣。

通過局部坐標(biāo)系下的虛位移方程得到局部坐標(biāo)系下對應(yīng)的結(jié)構(gòu)反力，最后通過逆變換得到對應(yīng)的整體坐標(biāo)系下的結(jié)構(gòu)反力：

FG=T-1R(TuG)

(10)

式中：FG為整體坐標(biāo)系下的節(jié)點(diǎn)反力，R的表達(dá)式見式(6)。

1.2 時(shí)間離散

傳統(tǒng)的平均加速度方法的時(shí)間離散采用梯形法則：

(11)

式中：u為節(jié)點(diǎn)位移，v為節(jié)點(diǎn)速度，a為節(jié)點(diǎn)加速度，由于非線性體系下平均加速度方法無法滿足能量守恒，因此需要其進(jìn)行修正。引入?yún)?shù)β對速度-加速度關(guān)系進(jìn)行修正，思路如下：

vi+1=vi+Δtaβ

(12)

修正項(xiàng)定義如下：

(13)

ri+β=R[βxi+(1-β)xi+1]

ri+1-β=R[(1-β)xi+βxi+1]

(14)

式中：f為外荷載。

假定計(jì)算步長內(nèi)的外荷載為常數(shù)，將其視為有勢力并考慮為系統(tǒng)的勢能的一部分，并定義系統(tǒng)能量增量為：

(15)

在無阻尼彈性系統(tǒng)中能量增量為零，將(13)和(11)第一式代入式(15)中得到文獻(xiàn)[17]給出的能量約束方程：

(16)

文獻(xiàn)[17]給出的變形能表達(dá)式在實(shí)際計(jì)算中運(yùn)用比較困難，本文通過變形能密度進(jìn)行體積分的方式計(jì)算變形能，并證明了兩種計(jì)算表達(dá)式等價(jià)，證明如下：

(17)

式中P為單元彈性勢能，根據(jù)虛位移原理：

(18)

代入式(16)，得到變形能表達(dá)式

(19)

因此單元的彈性勢能可以通過式(8)計(jì)算得到，整體的彈性勢能的計(jì)算式如下：

(20)

式中：Ept為結(jié)構(gòu)整體彈性勢能，函數(shù)P的表達(dá)式見式(8)，n為單元數(shù)目。結(jié)合式(11)、式(15)以及平衡方程得到基于能量一致的無條件穩(wěn)定逐步積分方法：

(21)

在考慮阻尼的情況下真實(shí)系統(tǒng)存在耗散，自身并非能量守恒。本文提出的算法可以實(shí)現(xiàn)能量一致。能量一致的定義為：通過數(shù)值算法得到的系統(tǒng)能量能夠與真實(shí)系統(tǒng)一樣具有能量耗散的能力，則稱該算法具有能量一致的特性。

本文基于能量一致的概念，給出考慮阻尼的能量一致方法如下：

(22)

式中c為阻尼矩陣,參數(shù)aβ為：

(23)

(24)

將其代入能量增量定義式(15)中得到數(shù)值計(jì)算的系統(tǒng)能量增量：

(25)

按照式(15)給出的能量增量的定義，黏滯阻尼系統(tǒng)的實(shí)際能量增量為

(26)

對比式(25)可以看出，本文提出的算法在考慮阻尼的情況下可以實(shí)現(xiàn)能量一致，即算法能夠保證系統(tǒng)能量與真實(shí)系統(tǒng)一樣具有能量耗散特性。

2 數(shù)值算例

本部分通過三個(gè)算例驗(yàn)證算法的性能，計(jì)算程序平臺(tái)為MATLAB。

對于非線性體統(tǒng)，系統(tǒng)能量的失穩(wěn)是判斷其算法失穩(wěn)的標(biāo)準(zhǔn)之一。對于有阻尼系統(tǒng)，由于阻尼項(xiàng)在系統(tǒng)中起到耗能的作用，因此在算例中不容易觀察到計(jì)算失穩(wěn)的現(xiàn)象，因此本文給出的算例都是無阻尼系統(tǒng)的算例。

本文不采用傳統(tǒng)的牛頓迭代方法求解時(shí)間離散后的非線性方程組(21)：

g(vi+1,β)=0

(27)

因?yàn)閷τ诰€彈性體系：

ri+β+ri+1-β=k(xi+1+xi)

(28)

在第i時(shí)刻位移、速度、加速度已知的情況下，i+1時(shí)刻的位移、速度、加速度可以通過式(21)的前三項(xiàng)確定，即先給定一個(gè)迭代初始值β，然后通過式(21)的前三項(xiàng)得到對應(yīng)i+1時(shí)刻的值：

(29)

計(jì)算得到的值代入式(21)第四式中，判斷計(jì)算計(jì)算結(jié)果是否滿足能量守恒，如果滿足則結(jié)束方程組求解，如果不滿足則將式(29)計(jì)算得到的結(jié)果代入到式(21)中，計(jì)算出一個(gè)新的迭代值β：

β=B(xi,vi,vi+1,xi+1)

(30)

通過重復(fù)上述計(jì)算即可得到非線性方程組(21)的一組解。需要注意的是，參數(shù)β可能存在多解的情況(在線性系統(tǒng)中為任意解)，但是在計(jì)算過程中只需要找到滿足非線性方程的一組解即可。

2.1 懸臂梁靜力分析

懸臂梁長l=1 m，寬b=0.1 m，高h(yuǎn)=0.1 m，彈性模量E=206 GPa，自由端施加垂直桿長的集中荷載F。沿桿長方向的積分采用高斯積分，取5個(gè)積分點(diǎn)。采用牛頓迭代求解非線性方程組，迭代誤差限為10-6m。無量綱荷載Fl/EI分別取為0.2、0.4、0.6、0.8、1.0、1.2、1.4、1.6、1.8、2.0，其中EI為抗彎剛度。按本文方法得到的端部撓度與陳至達(dá)[20]通過橢圓積分得到的精確解見圖2，可以觀察到當(dāng)單元數(shù)目增加的時(shí)候，數(shù)值解會(huì)趨近于精確解，由此證明了本文提出的有限元模型的合理性。從圖2還可以看出，計(jì)算誤差不僅與單元數(shù)目有關(guān)，也與梁的變形大小有關(guān)。當(dāng)單元數(shù)量為10個(gè)時(shí)，部分不同荷載(或位移)下的精度見表1。

另一方面，桿在大變形下的主要特征之一是其沿桿長方向會(huì)產(chǎn)生不均勻的軸向應(yīng)變。通過對上述懸臂梁結(jié)構(gòu)采用一個(gè)單元進(jìn)行靜力分析，得到了荷載作用下軸向應(yīng)變沿桿長方向分布曲線如圖3所示，可以看出當(dāng)荷載增大的時(shí)候應(yīng)變的不均勻性更顯著。

圖2 梁自由端豎向荷載-位移曲線

圖3 應(yīng)變沿桿長的分布

Fl/EI十個(gè)單元/m陳志達(dá)/m誤差/%0.20.066 10.066 30.40.40.1290.1311.60.60.1860.1923.20.80.2370.2494.81.00.2820.3026.72.00.4350.49312.8

2.2 懸臂梁自由振動(dòng)

懸臂梁采用2個(gè)單元構(gòu)造，每個(gè)單元長l=1 m，寬b=0.1 m，高h(yuǎn)=0.05 m，彈性模量E=100 GPa，采用集中質(zhì)量矩陣，質(zhì)量104kg，轉(zhuǎn)動(dòng)慣量103kg/m-2。

給予自由端豎向初始速度1 m/s。計(jì)算時(shí)長5.5 s，時(shí)間步長0.05 s。使用牛頓迭代求解非線性方程組，雅可比矩陣采用割線近似切線的方法進(jìn)行近似求解。能量方法初始β取為0.2，采用能量增量百分比控制，收斂準(zhǔn)則為總能量的萬分之一。計(jì)算了時(shí)間步長從0.05 s、0.06 s～0.15 s等十一個(gè)工況，能量方法在所有這些工況的計(jì)算結(jié)果均保持穩(wěn)定。當(dāng)時(shí)間增加到步長0.09 s的時(shí)候，平均加速度方法出現(xiàn)失穩(wěn)，步長為0.11 s的工況見圖6。圖7表明本文提出的方法計(jì)算結(jié)果穩(wěn)定并且很好地保證了系統(tǒng)的能量守恒。圖8表明參數(shù)β在遞推過程中大多數(shù)情況下為常數(shù)，對于該現(xiàn)象的解釋有待后續(xù)研究。

圖4 梁自由振動(dòng)示意圖

圖5 Δt=0.05 s梁自由端豎向位移時(shí)程曲線

圖6 Δt=0.11 s梁自由端豎向位移時(shí)程曲線

當(dāng)計(jì)算步長為0.05 s、總持時(shí)為5.5 s時(shí)，平均加速度方法耗時(shí)0.454 s，共迭代552次；能量方法耗時(shí)0.524 s，共迭代713次。能量方法耗時(shí)稍長。但是，當(dāng)能量方法步長取0.11 s時(shí)，則耗時(shí)降至0.331 s，迭代次數(shù)降至408次。因此，可以通過增加計(jì)算步長來減少能量方法計(jì)算代價(jià)。電腦配置為8.00 GB內(nèi)存，處理器為Inter(R)Core(TM)i7-7700HQ CPU@ 2.80 GHz。

圖7 Δt=0.11 s總能量時(shí)程曲線

圖8 Δt=0.11 s參數(shù)β-時(shí)間圖

圖9 剛架受力示意圖

2.3 剛架強(qiáng)迫振動(dòng)

剛架的梁柱單元幾何尺寸相同，每個(gè)構(gòu)件按一個(gè)單元計(jì)算，每個(gè)單元長l=1 m，寬b=0.1 m，高h(yuǎn)=0.05 m，彈性模量E=100 GPa， 采用集中質(zhì)量矩陣，質(zhì)量104kg，轉(zhuǎn)動(dòng)慣量103kg/m-2。剛架左端節(jié)點(diǎn)施加幅值為2 000 kN、頻率為2.5 Hz的水平向正弦荷載。當(dāng)時(shí)間步長為0.04 s的時(shí)候觀察到平均加速度方法失穩(wěn)，見圖11。能量方法通過能量守恒保持了其穩(wěn)定性，而平均加速度方法的能量波動(dòng)較大，見圖12。

圖10 Δt=0.04 s剛架1節(jié)點(diǎn)轉(zhuǎn)角時(shí)程曲線

圖11 Δt=0.04 s系統(tǒng)總能量時(shí)程曲線

3 結(jié) 論

本文通過獨(dú)立插值的方式，考慮了大轉(zhuǎn)角并用工程應(yīng)變描述軸向應(yīng)變，提出了平面大變形歐拉梁的有限元模型。從能量守恒的角度出發(fā)，采用單參數(shù)修正方法，實(shí)現(xiàn)了歐拉梁單元的無條件穩(wěn)定逐步積分積分算法。算例驗(yàn)證了本文方法與平均加速度法相比具有更優(yōu)的數(shù)值穩(wěn)定性，體現(xiàn)了能量方法在大變形歐拉梁動(dòng)力分析上的優(yōu)異性能。

目前關(guān)于修正參數(shù)的研究還不完善，主要表現(xiàn)在兩個(gè)方面：

(1)非線性方程組求解方法依賴于迭代初值的選擇，而參數(shù)的迭代初值的選擇有待進(jìn)一步研究。

(2)通過算例2可以看出，在逐步積分過程中修正參數(shù)趨于一個(gè)常數(shù)，如果能利用該現(xiàn)象對現(xiàn)有的非線性方程求解方法進(jìn)行改良，則計(jì)算代價(jià)會(huì)將進(jìn)一步減少。