基于能量守恒的平面大變形歐拉梁逐步積分算法

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

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

逐步積分算法是求解結構動力響應的重要手段之一，它同時也被運用到混合試驗這類試驗與數值計算相結合的試驗中[1-2]。而無條件穩定的逐步積分算法能夠保證混合試驗這類極度依賴算法穩定性的試驗順利進行，因此在此類依賴算法穩定性的計算中具有巨大優勢。

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

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

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

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

1 能量方法推導

1.1 空間離散

假定梁的初始構型和當前構型，如圖1所示。

引入一個無量綱參數s∈[0,1]，它表示每個點對應在梁上的位置，s(0)代表梁左端點，s(1)代表梁右端點，令初始構型梁上點的坐標與參數s滿足函數關系：

圖1 初始構型與當前構型圖

(1)

式中：下標o表示初始構型，n表示當前構型，l為單元初始長度，梁單元的位移場滿足式(2)所示函數關系：

(2)

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

(3)

由于在運用虛位移方程的時候，本文采用基于初始構型的虛位移方程，因此曲率的表達式為：

(4)

軸向應變的表達式為：

(5)

通過對初始構型使用虛位移原理得到結構反力表達式：

(6)

對于線彈性本構，其結構反力表達式為：

(7)

式中：E為單元彈性模量，I為截面慣性矩，A為單元截面面積。由此得到梁單元的等效節點反力與節點位移之間的函數關系。同時為了運用能量方法，需要給出計算單元的變形能的計算表達式。對于矩形截面，變形能EP與節點位移的函數關系為：

(8)

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

上述推導是在局部坐標進行的，對于多單元的整體坐標，可以通過如下轉換得到整體坐標下的位移對應的結構反力。整體坐標下單元節點位移uG與局部坐標單元節點位移uL存在轉換關系：

uL=TuG

(9)

式中：T為坐標轉換矩陣。

通過局部坐標系下的虛位移方程得到局部坐標系下對應的結構反力，最后通過逆變換得到對應的整體坐標系下的結構反力：

FG=T-1R(TuG)

(10)

式中：FG為整體坐標系下的節點反力，R的表達式見式(6)。

1.2 時間離散

傳統的平均加速度方法的時間離散采用梯形法則：

(11)

式中：u為節點位移，v為節點速度，a為節點加速度，由于非線性體系下平均加速度方法無法滿足能量守恒，因此需要其進行修正。引入參數β對速度-加速度關系進行修正，思路如下：

vi+1=vi+Δtaβ

(12)

修正項定義如下：

(13)

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

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

(14)

式中：f為外荷載。

假定計算步長內的外荷載為常數，將其視為有勢力并考慮為系統的勢能的一部分，并定義系統能量增量為：

(15)

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

(16)

文獻[17]給出的變形能表達式在實際計算中運用比較困難，本文通過變形能密度進行體積分的方式計算變形能，并證明了兩種計算表達式等價，證明如下：

(17)

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

(18)

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

(19)

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

(20)

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

(21)

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

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

(22)

式中c為阻尼矩陣,參數aβ為：

(23)

(24)

將其代入能量增量定義式(15)中得到數值計算的系統能量增量：

(25)

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

(26)

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

2 數值算例

本部分通過三個算例驗證算法的性能，計算程序平臺為MATLAB。

對于非線性體統，系統能量的失穩是判斷其算法失穩的標準之一。對于有阻尼系統，由于阻尼項在系統中起到耗能的作用，因此在算例中不容易觀察到計算失穩的現象，因此本文給出的算例都是無阻尼系統的算例。

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

g(vi+1,β)=0

(27)

因為對于線彈性體系：

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

(28)

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

(29)

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

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

(30)

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

2.1 懸臂梁靜力分析

懸臂梁長l=1 m，寬b=0.1 m，高h=0.1 m，彈性模量E=206 GPa，自由端施加垂直桿長的集中荷載F。沿桿長方向的積分采用高斯積分，取5個積分點。采用牛頓迭代求解非線性方程組，迭代誤差限為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為抗彎剛度。按本文方法得到的端部撓度與陳至達[20]通過橢圓積分得到的精確解見圖2，可以觀察到當單元數目增加的時候，數值解會趨近于精確解，由此證明了本文提出的有限元模型的合理性。從圖2還可以看出，計算誤差不僅與單元數目有關，也與梁的變形大小有關。當單元數量為10個時，部分不同荷載(或位移)下的精度見表1。

另一方面，桿在大變形下的主要特征之一是其沿桿長方向會產生不均勻的軸向應變。通過對上述懸臂梁結構采用一個單元進行靜力分析，得到了荷載作用下軸向應變沿桿長方向分布曲線如圖3所示，可以看出當荷載增大的時候應變的不均勻性更顯著。

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

圖3 應變沿桿長的分布

Fl/EI十個單元/m陳志達/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 懸臂梁自由振動

懸臂梁采用2個單元構造，每個單元長l=1 m，寬b=0.1 m，高h=0.05 m，彈性模量E=100 GPa，采用集中質量矩陣，質量104kg，轉動慣量103kg/m-2。

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

圖4 梁自由振動示意圖

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

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

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

圖7 Δt=0.11 s總能量時程曲線

圖8 Δt=0.11 s參數β-時間圖

圖9 剛架受力示意圖

2.3 剛架強迫振動

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

圖10 Δt=0.04 s剛架1節點轉角時程曲線

圖11 Δt=0.04 s系統總能量時程曲線

3 結 論

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

目前關于修正參數的研究還不完善，主要表現在兩個方面：

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

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