一種基于應變插值大變形梁單元的剛-柔耦合動力學分析

張志剛， 齊朝暉, 吳志剛,2

(1.大連理工大學工業裝備結構分析國家重點實驗室， 遼寧 大連 116024；2.大連理工大學航空航天學院， 遼寧 大連 116024)

一種基于應變插值大變形梁單元的剛-柔耦合動力學分析

張志剛1， 齊朝暉1, 吳志剛1,2

(1.大連理工大學工業裝備結構分析國家重點實驗室， 遼寧 大連 116024；2.大連理工大學航空航天學院， 遼寧 大連 116024)

柔性多體系統動力學中的“動力剛化”現象起因于變形間的耦合。一次近似模型成功地解決了小變形情況下的剛柔耦合建模問題，但在大變形情況下則需要考慮更多的耦合效應。選取表征梁彎曲應變的曲率和軸向應變作為單元參數進行離散；在大變形大轉動基礎上得到了單元兩端節點運動學參數的遞推關系，構造出了能夠自動計及“動力剛化項”且適用于大變形剛柔耦合動力學分析的平面梁單元。最后采用所提應變插值單元求解了包含大變形和剛柔耦合動力學柔性梁的數值算例，驗證了算法的正確性和有效性。

剛柔耦合； 梁單元； 應變插值； 動力剛化

引 言

梁的剛柔耦合動力學分析是近年來力學、機械及航空航天領域研究的一個熱點和難點[1-12]。柔性梁的運動包含大范圍剛體運動與變形運動的耦合，且描述梁變形所需的形心位移和截面角之間也存在耦合關系。這些耦合的作用機理非常復雜，給做大范圍運動梁的剛柔耦合動力學的建模和求解帶來了困難。

傳統的柔性多體系統動力學混合坐標法直接套用了結構力學小變形假設，Kane[1]在用其對旋轉柔性梁進行求解時發現高速轉動時數值結果發散，表現為隨著轉速增加梁變得更加柔軟，這與事實相違背，并首次提出了“動力剛化”的概念。這一現象引發了國內外廣大學者的關注，此后大量研究圍繞剛體運動與彈性運動的耦合及“動力剛化”問題展開：洪嘉振[2-6]指出，傳統混合坐標法在小變形基礎上忽略了縱向位移對軸向位移的影響是一種零次近似模型，低速轉動忽略掉變形耦合對仿真結果影響不大，但在高速轉動情況下這種近似將會產生完全錯誤的結果；其提出小變形情況下在軸向位移表達式中補充縱向位移二階項的所謂一次近似模型，用該模型進行了中心剛體旋轉柔性梁及末端帶集中質量的旋轉柔性梁動力學特性，并進一步通過實驗驗證了一次近似模型能夠捕捉“動力剛化項”。章定國[7]在一次近似模型基礎上對平面運動柔性梁進行了進一步研究，指出大范圍運動與彈性運動的耦合不僅會產生“動力剛化”現象，還會產生“動力柔化”現象，并對旋轉內接懸臂梁的“動力柔化”效應進行了分析。和興鎖等[8]在Shi P等[9]研究基礎上提出一次近似模型未能充分考慮幾何非線性變形對橫向變形的影響，因此在縱向與橫向位移中增加了新的變形二次耦合項。

上述研究大都局限于單元本身為小變形情況，針對包含大變形的柔性梁剛柔耦合問題，劉錦陽[10]采用絕對坐標法對旋轉柔性梁進行了研究，并與一次近似模型進行了對比，發現隨著梁剛度降低一次近似模型結果發散，而絕對坐標梁單元能夠得到收斂的結果。陳思佳[11]在一次近似模型基礎上采用模態坐標對梁進行離散，通過保留了動力學中方程中的變形耦合量相關的高階項，對考慮大變形的平面梁剛柔耦合動力學分析做了初步研究，但是通過模態坐標離散的求解精度和效率取決于模態集的選取方式。

現有研究大都選取節點位移作為有限元離散參量，在小變形情況下根據伯努利梁變形耦合關系添加位移高階耦合項來構造了一次近似模型或者高次近似模型。而直接構造滿足大變形大轉動的變形耦合位移梁單元仍然十分困難和復雜。本文在對曲率插值梁單元研究[12]的基礎上，選取梁的彎曲應變(曲率)和軸向應變作為單元參數進行離散，從而得到了自動計及位移耦合關系的梁單元。采用本單元對平面梁大范圍運動剛柔耦合問題進行了研究，所得結果與相關工作[13-16]吻合，從而驗證了所構造單元的正確性和有效性。

1 大范圍運動梁單元虛功率方程

1.1 梁單元運動學描述

平面梁單元如圖1所示，整體慣性坐標系為{g1g2}，在梁左端截面建立固結在截面上的單元坐標系{e1e2}，第一根軸e1指向梁長方向，初始時刻與整體坐標系第一根軸夾角φ0。在每個梁截面建立隨截面轉動的截面坐標系{t1t2}，變形前截面坐標系與單元坐標系方向保持一致。

圖1 大范圍運動平面梁單元Fig.1 A plane beam element undergoes large range of motion

如圖1所示，弧長s處截面形心矢徑為

(1)

式中r1為單元坐標系原點相對整體坐標系原點矢徑，ρ為弧長s處截面形心相對于單元坐標系原點矢徑。取出原始弧長為ds的一段梁微元，變形后幾何關系

(2)

(3)

其中 (·)′=d(·)/ds，θ為截面相對于單元坐標系的轉角。

單元左端截面相對于整體坐標系的轉角φ1表征了單元大范圍剛體轉動，截面相對于單元坐標系的轉角θ代表了單元的彈性變形，因此單元域內截面相對于整體坐標系的轉角為

(4)

1.2 梁單元應變場的離散

對于不超過彈性范圍的平面梁，其本構關系可用截面坐標系中力的分量描述為

(5)

式中F為截面軸向力；M為彎矩；EA,EIz分別為截面拉壓剛度和彎曲剛度。梁的應變虛功率

(6)

由上式可知，梁單元應變虛功率只與曲率和軸向應變有關，而不受單元位移和轉動大小的影響，因此選取梁曲率和軸向應變進行離散可以方便地得到單元應變虛功率和簡潔的單元剛度矩陣。

(7)

其中

(8)

其中ξ=s/l，梁曲率速率虛變分由式(7)可得

(9)

將式(7)，(9)及軸向應變ε是常量的假設代入單元應變虛功率式(6)

(10)

式中K為單元剛度陣

(11)

1.3 單元內運動學遞推

單元域內截面的轉角可由曲率插值函數式(7)積分得到

(12)

其中

(13)

轉角θ為截面相對于單元轉角，它完全由單元的彎曲變形確定，將其代入式(12)可得截面相對整體坐標系轉角為

(14)

由梁截面形心坐標弧長導數式(3)，對弧長s積分并注意到軸向應變ε為常量，可以得到截面形心坐標

(15)

從這兩個式子可以看出截面形心坐標與轉角和軸向應變高度耦合。需要說明的是傳統混合坐標零次近似模型及一次近似模型都假設轉角為小量，即：θ?0，從而將三角函數近似為：sinθ≈θ，cosθ≈1-θ2/2，并由式(3)得

(16)

將其帶入式(15)得

(17)

上式推導過程中忽略掉了轉角θ的二次以上項及θ和軸向應變ε的耦合項，式中加下劃線部分為一次近似模型在零次模型基礎上增加的橫向位移對軸線變形的二次耦合項。由此可見混合坐標一次近似模型不適用于包含大轉動的大變形情況。

本文不作小轉角假設，基于此推導的單元完整地保留了變形間的耦合關系，不僅適用于小轉角情況，同樣也能滿足大轉角情況。由式(12)可知截面轉角θ為弧坐標s的二次式，因此式(15)兩個積分不存在有理形式解，本文由高斯積分得到數值解。

直接由式(15)就可得到截面形心坐標場，但是在推導質量陣和廣義力時將會出現二重數值積分，從而會引起公式推導和數值求解上的不便。再次考察單元坐標系下截面形心坐標x,w在端部滿足

(18)

(19)

式中

(20)

(21)

弧長s處截面形心坐標為

(22)

其中單元坐標系與整體坐標系轉換矩陣

(23)

式中φ0為單元初始方位角。

對式(22)求時間導數得到截面形心速度

(24)

(25)

根據單元坐標系下截面形心坐標式(19)，可得

(26)

其中

(27)

(28)

將式(26)代入得到截面形心速度

(29)

對速度取變分

(30)

將式(26)代入式(25)得到截面形心加速度

(31)

其中加速度余量

(32)

對截面轉角式(14)求導可得截面轉動角速度及虛角速度

(33)

角加速度

(34)

為了方便推導，本文將整體坐標系坐標下截面形心坐標r和轉角φ組成列陣為q，兩端面運動學參數

(35)

由式(14)，(29)，(31)，(33)，(34)，單元兩端截面運動學遞推關系可以用矩陣表示為

(36)

其中

(37)

(38)

1.4 梁單元慣性力虛功率

根據伯努利梁理論，平面梁可以看作由無數做平面運動的剛性截面組成，其慣性力虛功率為

(39)

其中ρ代表梁的線密度，Iz代表梁截面慣性矩。

將式(30)～(34)代入上式

(40)

質量矩陣和廣義力為

(41)

(42)

其中

(43)

(44)

質量陣和廣義力中常系數矩陣

(45)

1.5 梁單元外力虛功率

梁單元上作用的外力虛功率為

(46)

式中p(s)為分布載荷；fi,mi分別為單元節點處所受集中力和力矩。將式(30),(33)代入上式

(47)

其中

(48)

(49)

對于考慮重力的情況，可將重力看做一種特殊的分布力，其對應廣義力為

(50)

(51)

另外，還可以將結構阻尼當做外力處理，本文采用文獻[13]阻尼模型，單元結構阻尼力虛功率為

(52)

式中η為阻尼系數。對比單元應變虛功率式(10)，可將機構阻尼力虛功率表示為

(53)

其中廣義結構阻尼力列陣為

(54)

2 平面梁剛柔耦合動力學方程

2.1 單元間運動學遞推

(1)單元坐標系原點選在靠近起始點一端節點，第一根軸方向背離起始節點指向第二個節點;

(2)沿背離起始節點方向單元編號遞增。

單元連接拓撲關系通過定義鄰接單元數組L描述：L為1×n行向量，其中L(i)代表與第i號單元相鄰并靠近起始節點1的那一側單元編號；規定與起始節點相連單元的鄰接單元為0。

如圖2所示，下面單元標號滿足規則標號，箭頭標出了單元坐標系原點所在節點和第一根坐標軸方向。

圖2 梁單元規則標號Fig.2 Regularly labeling of beam elements

按照上述約定，圖示由三個單元組成的梁的鄰接單元數組為：L=[0 1 2]。

利用鄰接單元數組L，根據單元內運動學遞推關系式(36)，可將i號單元坐標系所在節點運動參量與起始節點關系表示為

(55)

其中

(56)

(57)

2.2 梁系統剛柔耦合動力學方程

彈性體虛功率方程表述為：彈性體所受外力虛功率等于慣性力虛功率和應變虛功率之和。劃分為n個單元柔性梁的虛功率方程為

(58)

將梁單元應變虛功率式(10)、慣性力虛功率式(40)、外力虛功率式(47)及單元間運動學遞推關系式(55)代入上式得到

其中質量矩陣

(60)

廣義力列陣

(61)

由梁系統虛功率方程式(59)中虛變分獨立，其系數為零可得平面梁動力學方程

(62)

3 數值算例

在這一節中采用本文單元編寫了動力學分析程序，分別選取了懸臂梁純彎曲、大范圍運動已知的旋轉柔性梁及包含大變形的自由下落柔性單擺3個算例進行了求解，并將所得結果與相關文獻[10-11,14-15]進行了分析比較。

3.1 懸臂梁受梁端集中彎矩作用

懸臂梁端受集中彎矩M作用，梁長度L=10 m，梁端線位移u，v如圖3所示。

圖3 懸臂梁受集中彎矩作用Fig.3 Cantilever beam undergoes bending load

由于存在結構阻尼，梁最終會處于靜平衡狀態，這里采用本文算法將整個梁劃分為5個單元，為了減少懸臂梁達到平衡狀態所用時間，這里選取阻尼比η=0.1進行仿真。在不同彎矩作用下梁最終的變形狀態如圖4所示。

圖4 不同彎矩作用下懸臂梁的大變形Fig.4 Large deformation of the cantilever under bending moments

下面給出懸臂梁在彎矩M=2πEIz/L作用下梁末端位移隨時間變化曲線如圖5所示。

圖5 懸臂梁末端受集中彎矩的位移Fig.5 The tip displacements of the cantilever with an end bending load

從位移時間曲線圖可以看出，在較大結構阻尼作用下懸臂梁經過兩個周期就達到了最終的平衡狀態，且僅用5個單元就得到了與解析解高度吻合的結果，由此表明本文單元在處理包含大變形問題的優越性。

3.2 旋轉柔性梁大范圍運動已知

大范圍運動已知的旋轉柔性梁已被很多學者[14-15]用來考察所提方法補充“動力剛化項”的能力。如圖6所示，柔性梁物理參數選取與文獻[14]相同：梁長L=10 m截面面積A=1 m2，截面慣性矩Iz=5×10-4m4，彈性模量E=2.8×107Pa，密度ρ=1.2 kg/m3。

圖6 旋轉柔性梁Fig.6 A flexible beam on a rotating base

柔性梁繞轉動中心的轉角φ(t)規律為

(63)

這里忽略梁結構阻尼，將梁劃分為5個單元，仿真時間為30s，得到梁自由端位移響應曲線如圖7所示。

圖7 自由端位移曲線Fig.7 Displacement curves at the free end

從圖7可以看出：梁端部位移隨著轉動速度增加而增大，當轉速增加到一定數值維持不變，軸向位移逐漸回落并穩定在5.14×10-4m，橫向位移在0附近小幅度波動。采用本文單元的仿真結果與文獻[14-15]完全相符，可以驗證算法對動力剛化問題的適用性。

3.3 自由下落柔性單擺

最后一個算例研究柔性單擺在重力作用下自由下落問題，單擺初始時刻水平放置如圖8所示。

圖8 柔性單擺Fig.8 The flexible pendulum

這個算例的大范圍剛體運動是未知的，是個典型的大范圍運動與變形運動耦合的問題。針對這個算例劉錦陽[10]和陳思佳[11]分別采用不同建模方法進行仿真求解，并與混合坐標法一次近似模型進行了對比。通過選取較低的材料彈性模量使得彈性變形不再是小量，以此來驗證方法對大變形問題的適用性。柔性擺的長度L=1.8 m，E=6.895×109Pa，A=2.5 cm2，Iz=0.13 cm4，密度ρ=2 766.67 kg/m3。

這個算例中將單擺劃分為5個單元，以本文方法和絕對節點坐標方法[10]在Matlab環境中選用剛性微分方程組求解器ode23tb，數值積分參數設定為：絕對誤差限AbsT=0.1，相對誤差限RelT=0.01，在同一臺微機上進行仿真2.5 s，所用CPU時間如表1所示。

表1 絕對坐標方法與本文方法CPU時間對比

Tab.1 The CPU time comparison between absolute coordinate method and the proposed method

方法文獻[10]絕對坐標法本文方法CPU時間225．05s6．29s

兩種方法所得單擺末端位移時間曲線如圖9所示。

圖9 柔性單擺的末端位移Fig.9 The tip deformation of the flexible pendulum

根據文獻[10-11]結果，混合坐標法一次近似模型仿真結果發散，表明其不適用與包含大變形的梁的剛柔耦合動力學分析。而本文采用應變插值所得結果與文獻[10-11]相吻合，具有較好的計算精度。

從兩種方法所用CPU時間對比可以看出，本文單元求解效率優于文獻[10]絕對節點坐標單元，也具有很好的計算效率。

4 結 論

本文通過直接選取單元應變進行插值構造了一種不受單元變形大小限制，適用于剛柔耦合動力學分析的平面梁單元。梁變形的位移及轉角之間的耦合關系通過幾何方程積分得到，并由此給出了節點運動學遞推關系。本文所提梁單元有如下特點：

(1)選取單元應變進行離散不受單元剛體運動的影響，得到了簡潔的單元剛度矩陣及節點力；

(2)根據應變場由幾何方程積分得到梁截面形心坐標和轉角耦合關系，自動計及了“動力剛化項”；

(3)單元不受變形大小的限制，適用于包含大變形梁的剛柔耦合動力學分析。

[1] Kane T R, Ryan R, Banerjee A K. Dynamics of a cantilever beam attached to a moving base[J]. Journal of Guidance, Control, and Dynamics, 1987, 10(2): 139—151.

[2] 楊輝, 洪嘉振, 余征躍. 動力剛化問題的實驗研究[J]. 力學學報, 2004, 36(1): 118—124.

YANG Hui, HONG Jiazhen, YU Zhengyue. Experimental investigation on dynamic stiffening phenomenon[J]. Acta Mechanica Sinica, 2004，36(1)：119—124.

[3] 楊輝, 洪嘉振, 余征躍. 剛柔耦合建模理論的實驗驗證[J]. 力學學報, 2003, 35(2): 253—256.

YANG Hui, HONG Jia-zhen，YU Zheng-yue．Experiment validation of modeling theory for rigid-flexible coupling systems[J]．Acta Mechanica Sinica，2003，35(2)：253—256．

[4] Cai G P, Hong J Z, Yang S X. Dynamic analysis of a flexible hub-beam system with tip mass[J]. Mechanics Research Communications, 2005, 32(2): 173—190.

[5] Cai G P, Lim C W. Dynamics studies of a flexible hub-beam system with significant damping effect[J]. Journal of Sound and Vibration, 2008, 318(1): 1—17.

[6] 劉錦陽, 洪嘉振. 柔性梁的剛-柔耦合動力學特性研究[J]. 振動工程學報, 2002, 15(2): 194—198.

Liu Jinyang, Hong Jiazhen. Study on rigid-flexible coupling dynamic behaviour of flexible beam[J]. Journal of Vibration Engineering,2002, 15(2): 194—198.

[7] 方建士, 章定國. 旋轉內接懸臂梁的剛柔耦合動力學特性分析[J]. 物理學報, 2013，62(4): 305—311.

Fang Jianshi, Zhang Dingguo. Analyses of rigid-flexible coupling dynamic properties of a rotating internal cantilever beam[J]. Acta Physica Sinica,2013,62(4) : 305—311.

[8] 鄧峰巖, 和興鎖, 楊永鋒, 等. 計及非線性變形的剛-柔耦合動力學建模[J]. 機械強度, 2006, 28(6): 800—804.

Deng Fengyan, He Xingsuo, Yang Yongfeng, et al. Dynamics modeling for a rigid-flexible coupling system with onlinear deformation field[J]. Journal of Mechanical Strength， 2006, 28(6): 800—804

[9] Shi P, McPhee J, Heppler G R. A deformation field for Euler-Bernoulli beams with applications to flexible multibody dynamics[J]. Multibody System Dynamics, 2001, 5(1): 79—104.

[10]李彬, 劉錦陽. 大變形柔性梁系統的絕對坐標方法[J]. 上海交通大學學報, 2005, 39(5): 827—831.

Li Bin, Liu Jinyang. Application of absolute nodal coordination formulation in flexible beams with large deformation[J]. Journal of Shanghai Jiaotong University, 2005, 39(5):827—831．

[11]陳思佳, 章定國, 洪嘉振. 大變形旋轉柔性梁的一種高次剛柔耦合動力學模型[J]. 力學學報, 2013, 45(2): 251—256.

Chen Sijia, Zhang Dingguo, Hong Jiazhen. A high-order rigid-flexible coupling model of a rotating flexible beam under lagre deformation[J]. Acta Mechanica Sinica, 2013, 45(2): 251—256.

[12]張志剛, 齊朝暉, 吳志剛. 基于曲率插值的大變形梁單元[J]. 應用數學和力學, 2013 34(6): 620—629.

Zhang Zhigang, Qi Zhaohui, Wu Zhigang. A large deformation beam element based on curvature interpolation[J]. Applied Mathematics and Mechanics, 2013, 34(6): 620—629.

[13]Garcíd D, Valverde J, Dominguez J. An internal damping model for the absolute nodal coordinate formulation[J]. Nonlinear Dynamics, 2005, 42(4): 347—369.

[14]Simo J C, Vu-Quoc L. On the dynamics in space of rods undergoing large motions—a geometrically exact approach[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 66(2): 125—161.

[15]Zhao Z, Ren G. A quaternion-based formulation of Euler-Bernoulli beam without singularity[J]. Nonlinear Dynamics, 2012, 67(3):1 825—1 835.

Rigid-flexible dynamics analysis of a large deformation beam element based on interpolation of strains

ZHANGZhi-gang1,QIZhao-hui1,WUZhi-gang1,2

(1.State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China; 2.School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China)

The “dynamic stiffening” phenomenon in flexible multi-body system dynamics is due to the deformation coupling. The first-order approximation model has been successfully applied in the rigid-flexible coupling modeling of small deformation. However, it is found necessary to consider more deformation coupling effects in lager deformation cases. In this paper, the beam bending curvature and axial strain are selected as the element parameters. Then the recursion formulations are obtained for the kinematic parameters of two end nodes of an element based on theories of large deformation and finite rotation. A planar beam element used for large deformation rigid-flexible dynamics analysis is proposed, which can automatically take into account the “dynamic stiffening terms”. Finally, the validity and effectiveness of the proposed algorithm are verified through some numerical examples which involve the large deformations and rigid-flexible dynamics of beams.

rigid-flexible coupling; beam element; strain interpolation; dynamic stiffening

2013-12-26；

2014-06-12

國家自然科學基金資助項目(11372057)

O313.7

A

1004-4523(2015)03-0337-08

10.16385/j.cnki.issn.1004-4523.2015.03.001

張志剛(1984—),男，博士研究生。電話： 15326175369； E-mail:zhigangzhang@foxmail.com

齊朝暉(1964—),男，教授，博士生導師。E-mail: zhaohuiq@dlut.edu.cn