徐 圣，劉錦陽，余征躍

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

由多個(gè)柔性部件組成的多體系統(tǒng)稱為柔性多體系統(tǒng)，如大型空間展開機(jī)構(gòu)、風(fēng)力發(fā)電機(jī)葉片、直升飛機(jī)機(jī)翼等。在外載荷作用下，每個(gè)柔性部件可能會(huì)產(chǎn)生較大的彈性變形。梁被定義為某一方向上的尺寸比另外兩個(gè)方向大得多的構(gòu)件，目前對大變形空間柔性梁的建模和分析因?yàn)閹缀畏蔷€性的存在而變得具有挑戰(zhàn)性，在保證計(jì)算精度的前提下，如何提高數(shù)值計(jì)算的收斂性，實(shí)現(xiàn)大變形空間柔性梁系統(tǒng)高效的動(dòng)力學(xué)仿真是一個(gè)非常有意義的命題。

混合坐標(biāo)法［1］采用剛體大范圍運(yùn)動(dòng)的剛體坐標(biāo)與柔性體的節(jié)點(diǎn)坐標(biāo)(或模態(tài)坐標(biāo))為廣義坐標(biāo)，建立柔性多體系統(tǒng)剛－柔耦合的動(dòng)力學(xué)模型。傳統(tǒng)的混合坐標(biāo)法在建模過程中，直接套用了以線性小位移為基本前提的結(jié)構(gòu)動(dòng)力學(xué)假設(shè)，忽略了應(yīng)變和位移關(guān)系式中彈性變形的高次項(xiàng)，并采用模態(tài)縮減法提高計(jì)算效率。然而，在大變形情況下，線彈性理論和模態(tài)縮減法均失效，需要考慮幾何非線性項(xiàng)，混合坐標(biāo)法的優(yōu)勢就無法體現(xiàn)。

絕對節(jié)點(diǎn)坐標(biāo)法［2－4］是目前應(yīng)用較為廣泛的大變形柔性多體系統(tǒng)動(dòng)力學(xué)建模方法。該方法的特點(diǎn)是選取柔性梁各單元節(jié)點(diǎn)相對慣性基的位移和斜率作為廣義坐標(biāo)，建立動(dòng)力學(xué)方程。這樣就不需要引入大范圍運(yùn)動(dòng)變量和變量位移坐標(biāo)，廣義坐標(biāo)全部是全局坐標(biāo)，得到的質(zhì)量陣為常數(shù)陣。Omar等［5］進(jìn)一步考慮剪切效應(yīng)，建立二維梁的有限元模型。在平面梁的基礎(chǔ)上，Shabana等［6］將絕對節(jié)點(diǎn)坐標(biāo)法推廣到大變形的空間梁。但是，由于取每個(gè)節(jié)點(diǎn)的絕對位移坐標(biāo)和沿局部坐標(biāo)系三個(gè)方向的斜率坐標(biāo)為廣義坐標(biāo)，空間梁的每個(gè)節(jié)點(diǎn)有12個(gè)自由度，計(jì)算量較大，收斂速度緩慢。雖然考慮了橫向剪切和截面變形，但是彎曲問題的數(shù)值收斂性低于預(yù)期，存在著剪切鎖定［7－8］問題。

解決大變形問題的另一種方法是幾何精確建模方法［9］，該方法將空間中軸線的三個(gè)方向的位移和橫截面的三個(gè)方向的轉(zhuǎn)角作為廣義坐標(biāo)，建立動(dòng)力學(xué)方程，但是在對三個(gè)橫截面轉(zhuǎn)角進(jìn)行有限元離散時(shí)，如果用線性插值函數(shù)表示有限轉(zhuǎn)角，容易引起數(shù)值計(jì)算誤差。為了解決這個(gè)問題，需要對后一個(gè)節(jié)點(diǎn)相對前一個(gè)節(jié)點(diǎn)的轉(zhuǎn)角坐標(biāo)進(jìn)行線性插值，使建模過程復(fù)雜化。由于工程中常見的空間梁的厚度和寬度的比值很小，橫向剪切可以不計(jì)，中軸線的切線和橫截面法線基本垂直，只需要考慮拉伸變形，彎曲變形和扭轉(zhuǎn)變形，這樣，除了扭轉(zhuǎn)角，其余轉(zhuǎn)角均可用絕對位移關(guān)于局部坐標(biāo)的偏導(dǎo)數(shù)表示，從而避免了用線性插值函數(shù)對有限轉(zhuǎn)角進(jìn)行插值引起的精度問題。

和興鎖等［10］考慮梁的縱向、橫向、側(cè)向彎曲變形以及扭轉(zhuǎn)變形的耦合項(xiàng)，利用Lagrange方程建立具有大范圍運(yùn)動(dòng)和非線性變形的柔性梁的有限元?jiǎng)恿W(xué)方程。陳思佳等［11］在動(dòng)力學(xué)方程中保留與非線性耦合量相關(guān)的某些高階項(xiàng)，建立了做大范圍旋轉(zhuǎn)運(yùn)動(dòng)的中心剛體－柔性懸臂梁系統(tǒng)的高次剛?cè)狁詈蟿?dòng)力學(xué)模型。目前，現(xiàn)有大型工程軟件LS-DYNA，ABAQUS等均可以用于計(jì)算大變形問題，但是由于這些軟件在動(dòng)力學(xué)建模時(shí)，為了簡化公式，對方向余弦陣作了近似，需要取較多的單元才能保證收斂和計(jì)算精度。

為了準(zhǔn)確而高效地解決大變形空間梁的動(dòng)力學(xué)問題，本文基于幾何精確建模方法建立了大變形細(xì)長空間梁的幾何非線性動(dòng)力學(xué)模型，并通過動(dòng)力學(xué)實(shí)驗(yàn)驗(yàn)證建模理論的正確性，通過與LS-DYNA的數(shù)值對比驗(yàn)證了本文方法的快速收斂性。首先用空間曲梁中線上任意一點(diǎn)的3個(gè)絕對位置坐標(biāo)和橫截面的3個(gè)姿態(tài)角描述橫截面的位置和姿態(tài)，詳細(xì)推導(dǎo)了應(yīng)變和曲率和位置、姿態(tài)坐標(biāo)的關(guān)系，在此基礎(chǔ)上基于中線切線與橫截面法線重合的假設(shè)，不計(jì)橫向剪切變形，將節(jié)點(diǎn)廣義坐標(biāo)數(shù)從6個(gè)減少為4個(gè)，簡化了動(dòng)力學(xué)模型，從而減小了數(shù)值仿真所需的計(jì)算時(shí)間。在大變形情況下，設(shè)計(jì)氣浮臺(tái)和柔性空間梁系統(tǒng)的剛－柔耦合動(dòng)力學(xué)實(shí)驗(yàn)平臺(tái)，通過對氣浮臺(tái)和梁系統(tǒng)的剛－柔耦合實(shí)驗(yàn)驗(yàn)證了本文的幾何非線性空間梁動(dòng)力學(xué)模型的準(zhǔn)確性。此外，就同一算例，分別用幾何非線性空間梁建模方法和現(xiàn)有大型工程軟件(LS-DYNA)進(jìn)行仿真計(jì)算，獲得一致的結(jié)果。進(jìn)一步將本文方法仿真計(jì)算所需的時(shí)間與LS-DYNA進(jìn)行比較，驗(yàn)證了本文方法具有良好的單元收斂性和高效性。本文方法的特點(diǎn)在于利用較少的梁單元即可獲得足夠高的計(jì)算精度，可有效減少計(jì)算耗時(shí)。

1 應(yīng)變和曲率的推導(dǎo)

作大范圍運(yùn)動(dòng)的空間柔性曲梁如圖1所示。建立三個(gè)相關(guān)的空間坐標(biāo)系:O－e1e2e3為固定的參考坐標(biāo)系，P－為固結(jié)于變形前過梁中軸線上P點(diǎn)的橫截面坐標(biāo)系，P'－為固結(jié)于變形后過梁中軸線上P'點(diǎn)的橫截面坐標(biāo)系?；【€坐標(biāo)s為曲梁根部至中軸線參考點(diǎn)的弧長的長度。變形前橫截面上任意一點(diǎn)的坐標(biāo)陣為:

圖1 變形前后的空間曲梁Fig.1 A spatial curved beam before and after deformation

不計(jì)橫截面的翹曲變形，并假定橫截面的中心在中軸線上，變形后橫截面上任意一點(diǎn)的坐標(biāo)陣為:

式中，y，z是橫截面上任意一點(diǎn)相對中軸線參考點(diǎn)的截面坐標(biāo)，r0，r分別是中軸線上P點(diǎn)變形前后在固定參考坐標(biāo)系下的坐標(biāo)陣，r0=［r01r02r03］T，r=［r1r2r3］T，A0，A 分別是變形前后橫截面坐標(biāo)系相對參考坐標(biāo)的方向余弦陣。

設(shè)橫截面上任意一點(diǎn)的坐標(biāo)陣對自然坐標(biāo)s求導(dǎo)用‘表示，根據(jù)幾何精確梁理論，梁上任意一點(diǎn)處的應(yīng)變列陣為:

式中，ε1表示軸向拉伸應(yīng)變，ε2，ε3表示橫向剪切變形。將式(1)和(2)代入式(3)，得到:

進(jìn)一步改寫得到:

為簡化表達(dá)，令

為提高幾何非線性空間梁的計(jì)算效率，考慮減小整個(gè)廣義坐標(biāo)陣的規(guī)模。對于細(xì)長梁來說，其橫向剪切應(yīng)變可以忽略。假定曲梁中軸線的切向與橫截面法向一致，則有:

中軸線上對應(yīng)點(diǎn)的應(yīng)變列陣可簡化為

利用關(guān)系式(12)得到:

由于2和3可以用r和r'表示，橫截面的位置和姿態(tài)可以用r，r'和1描述，從而縮減了節(jié)點(diǎn)自由度。將式(15)～(17)代入式(9)，得到:

2 柔性梁的有限元?jiǎng)恿W(xué)方程

下面采用2節(jié)點(diǎn)Hermite函數(shù)對每個(gè)梁單元進(jìn)行插值。設(shè)第i個(gè)梁單元坐標(biāo)陣qi，總體廣義坐標(biāo)陣q，Bi為單元坐標(biāo)陣qi與總體坐標(biāo)陣q之間的轉(zhuǎn)換陣，qi=Biq，則

將式(23)和(24)代入式(14)，彈性力的虛功率為:

設(shè)FP是作用于中軸線上P點(diǎn)的外力在O－e1e2e3上的坐標(biāo)陣，外力的虛功率為:

將式(37)和(38)代入(36)，外力偶的虛功率為:

其中，Φq為約束方程關(guān)于q的導(dǎo)數(shù)陣，λ為與約束力(偶)相關(guān)的拉格朗日列陣。方程(41)需要與運(yùn)動(dòng)學(xué)約束方程聯(lián)合求解。

3 數(shù)值仿真與討論

3.1 大變形懸臂梁動(dòng)力學(xué)響應(yīng)研究

不計(jì)重力，考慮一端懸臂的細(xì)長柔性直梁。A、B分別是梁的左端點(diǎn)和右端點(diǎn)，在梁右端點(diǎn)B處施加沿著z軸正向的外力偶M，懸臂梁做大范圍運(yùn)動(dòng)，如圖2所示。

圖2 力偶作用下的柔性懸臂梁Fig.2 A flexible cantilevered beam applied with tip torque

梁的長度為1 m，橫截面為邊長2 mm×2 mm的正方形。梁的材料參數(shù)為:ρ=2 700 kg/m3，E=7×1010Pa，μ =0.3。

外力偶M的大小隨時(shí)間的變化規(guī)律為:

在本文的幾何非線性建?；A(chǔ)上，用廣義－α方法(譜半徑取0.5)對動(dòng)力學(xué)方程進(jìn)行數(shù)值求解，同時(shí)利用LS-DYNA對該懸臂梁進(jìn)行動(dòng)力學(xué)仿真，得到梁的右端點(diǎn)y方向的位移和速度隨時(shí)間(1～10 s)的變化曲線，如圖3和圖4所示?？梢钥闯?，本文方法與LS-DYNA的計(jì)算結(jié)果吻合。

圖3 梁端點(diǎn)的y方向位移Fig.3 Tip displacement of the tip point in y direction

圖4 梁端點(diǎn)的y方向速度Fig.4 Tip velocity of the tip point in y direction

表1為本文方法和LS-DYNA的仿真耗時(shí)比較。可以看出，本文方法使用10個(gè)梁單元(最少收斂單元數(shù))，計(jì)算時(shí)間為42 s;LS-DYNA使用45個(gè)梁單元(最少收斂單元數(shù))，計(jì)算時(shí)間為200 s?？梢姡瑤缀畏蔷€性方法在保證小誤差的同時(shí)，由于10個(gè)梁單元就可以收斂，具有更高的計(jì)算效率。

表1 本文方法與LS-DYNA的仿真耗時(shí)對比Tab.1 Comparison of time cost of the present formulation and LS-DYNA

3.2 柔性梁的剛－柔耦合動(dòng)力學(xué)研究

如圖5所示，氣浮臺(tái)B1繞z軸作定軸轉(zhuǎn)動(dòng)，矩形柔性梁B2固結(jié)在氣浮臺(tái)上。xyz為固結(jié)于地面的慣性基，x1y1z1為氣浮臺(tái)連體基，x2y2z2為柔性梁浮動(dòng)基。梁的浮動(dòng)基的z2方向與氣浮臺(tái)的轉(zhuǎn)軸z1的夾角為a=π/4，重力沿著 z1的負(fù)向。梁長1.45 m，寬0.1 m，厚0.002 5 m，體密度 2 766 kg/m3，彈性模量 6.9 × 1010Pa。氣浮臺(tái)的轉(zhuǎn)軸與梁的內(nèi)端點(diǎn)的距離為0.35 m，氣浮臺(tái)的轉(zhuǎn)動(dòng)慣量為J=9.557 kg·m2，觀測點(diǎn)C和D到A點(diǎn)的距離分別為0.3 m和1.2 m 。

圖5 氣浮臺(tái)剛－柔耦合實(shí)驗(yàn)平臺(tái)Fig.5 Air-floating rigid-flexible coupling experiment test bed

以q—=［θqT］T為系統(tǒng)廣義坐標(biāo)，q為梁上各節(jié)點(diǎn)的廣義坐標(biāo)陣。由于 B2在A點(diǎn)與B1固結(jié)，且兩物體連體基中向量x1與向量x2平行，因此，y1和z1均垂直于x2。y1與y2之間的夾角為a=π/4，z1與y2之間的夾角為a+π/2=3π/4，則系統(tǒng)的運(yùn)動(dòng)學(xué)約束方程為:

其中，A(θ)為氣浮臺(tái)連體基x1y1z1關(guān)于慣性基xyz的方向余弦陣。

初始時(shí)刻，矩形梁無變形地?cái)R置在傾斜α=45°的平板上，釋放后柔性梁在重力作用下產(chǎn)生彈性變形，氣浮臺(tái)同時(shí)發(fā)生轉(zhuǎn)動(dòng)。

用本文幾何非線性空間梁理論對上述系統(tǒng)進(jìn)行建模仿真，并通過實(shí)驗(yàn)測試系統(tǒng)響應(yīng)。在仿真計(jì)算中，考慮了結(jié)構(gòu)阻尼力，阻尼陣C=βM(根據(jù)振幅衰減率得到阻尼系數(shù)β為0.2)。通過角速度傳感器測得氣浮臺(tái)的角速度以及觀測點(diǎn)C和D的上表面的縱向應(yīng)變隨時(shí)間的變化曲線，將測得的實(shí)驗(yàn)數(shù)據(jù)與仿真結(jié)果進(jìn)行對比，如圖6圖7和圖8所示。

圖6 轉(zhuǎn)臺(tái)的角速度Fig.6 Angular velocity of the hub

圖7 C點(diǎn)上表面的軸向應(yīng)變Fig.7 Axial strain of point C on the upper surface

圖8 D點(diǎn)上表面的軸向應(yīng)變Fig.8 Axial strain of point D on the upper surface

由上述曲線圖可以看出，幾何非線性方法計(jì)算所得的氣浮臺(tái)角速度，以及觀測點(diǎn)C和D上表面的軸向應(yīng)變與實(shí)驗(yàn)測量值基本吻合。實(shí)驗(yàn)值的頻率略微低于計(jì)算值，這是由于空間梁與氣浮臺(tái)的固結(jié)不是理想剛體約束，在實(shí)驗(yàn)中有微小變形。通過數(shù)值對比，驗(yàn)證了本文幾何非線性方法的準(zhǔn)確性。

4 結(jié)論

(1)與LS-DYNA相比，本文幾何非線性空間梁建模方法在保證誤差較的同時(shí)具有較高的計(jì)算效率。首先，該方法將每個(gè)節(jié)點(diǎn)的坐標(biāo)數(shù)降為4個(gè)，從而縮減了梁的整體廣義坐標(biāo)規(guī)模。其次，數(shù)值解收斂所需要的單元數(shù)少于LS-DYNA，有利于減小計(jì)算量。因此，幾何非線性空間梁建模方法在計(jì)算效率上具有優(yōu)越性。

(2)設(shè)計(jì)了氣浮臺(tái)－空間梁系統(tǒng)剛－柔耦合動(dòng)力學(xué)實(shí)驗(yàn)平臺(tái)，通過動(dòng)力學(xué)仿真結(jié)果與實(shí)驗(yàn)值的比較，驗(yàn)證了本文幾何非線性空間梁建模方法的準(zhǔn)確性。

［1］Linkins P W.Finite element appendage equations for hybrid coordinate dynamic analysis［J］.International Journal of Solids＆ Structures，1972，8(5):709－731.

［2］ Shabana A A，Schwertassek R.Equivalance of the floating frame of reference approach and finite element formulations［J］，International Journal of Non-Linear Mechanics，1998，33(3):417－432.

［3］ Shabana A A，Hussein H A，Escalona JL.Application of the absolute nodal coordinate formulation to large rotation and large deformation problems ［J］， ASME Journal of Mechanical Design，1998，120(2):188－195.

［4］Berzeri M，Shabana A A.Development of simple models for the elastic forces in the absolute nodal coordinate formulation［J］.Journal of Sound and Vibration，2000，235(4):539－565.

［5］Omar M Z，Shabana A A.A two-dimensional shear deformable beam for large rotation and deformation problems ［J］.Journal of Sound and Vibration，2001，243(3):565－576.

［6］Shabana A A，Yakoub R Y.Three dimensional absolute nodal coordinate formulation for beam elements:theory［J］.ASME Journal of Mechanical Design，2001，123(4):606－613.

［7］Gerstmayr J，Shabana A A.Analysis of thin beams and cables using the absolute nodal coordinate formulation ［J］.Nonlinear Dynamics，2006，45(1－2):109－130.

［8］ Schwab A L，Meijaard J P.Comparison of three-dimensional flexible beam elements for dynamic analysis:Finite element method and absolute nodal coordinate formulation ［J］，ASME Journal of Computational Nonlinear Dynamics，2005，5(1):1－10.

［9］Bauchau O A.Flexible multibody dynamics［M］.Dodrecht:Springer，2011.

［10］和興鎖，鄧峰巖，吳根勇，等.對于具有大范圍運(yùn)動(dòng)和非線性變形的柔性梁的有限元?jiǎng)恿W(xué)建模［J］.物理學(xué)報(bào)，2010，59(1):25－29.HE Xing-suo， DENG Feng-yan， WU Gen-yong， et al.Dynamic modeling of a flexible beam with large overall motion and nonlinear deformation using the finite element method［J］.Acta Physica Sinica，2010，59(1):25 －29.

［11］陳思佳，章定國，洪嘉振.大變形旋轉(zhuǎn)柔性梁的一種高次剛－柔耦合動(dòng)力學(xué)模型［J］.力學(xué)學(xué)報(bào)，2012，45(2):251－256.CHEN Si-jia，ZHANG Ding-guo，HONG Jia-zhen.A highorder rigid-flexible coupling model of a rotating flexible beam under large deformation［J］.Chinese Journal of Theoretical and Applied Mechanics，2012，45(2):251 －256.