考慮剛-柔-熱耦合的板結構多體系統(tǒng)的動力學建模*

張煒華 劉錦陽

（上海交通大學工程力學系，上海 200240）

考慮剛-柔-熱耦合的板結構多體系統(tǒng)的動力學建模*

張煒華 劉錦陽?

（上海交通大學工程力學系，上海 200240）

對熱載荷作用下中心剛體與大變形薄板多體系統(tǒng)的動力學建模問題進行研究.基于Kirchhoff假設，從格林應變和曲率與絕對位移的非線性關系式出發(fā)，推導了非線性廣義彈性力陣，用絕對節(jié)點坐標法建立了大變形矩形薄板的有限元離散的動力學變分方程.為了考慮剛體姿態(tài)運動、彈性變形和溫度變化的相互耦合作用，推導了熱流密度與絕對節(jié)點坐標之間的關系式.引入系統(tǒng)的運動學約束方程，建立了中心剛體-矩形板多體系統(tǒng)的考慮剛-柔-熱耦合的熱傳導方程和帶拉格朗日乘子的第一類拉格朗日動力學方程.為了有效地提高計算效率，將改進的中心差分法和廣義-α法相結合，求解熱傳導方程和動力學方程，差分后的方程通過牛頓迭代法耦合求解.對剛-柔耦合和剛-柔-熱三者耦合兩種模型的仿真結果進行比較表明，剛體運動對溫度梯度和熱變形的影響顯著.此外，本文建模方法考慮了幾何非線性項，因此也考慮了熱膨脹引起的軸向變形對橫向變形的影響.

大變形薄板多體系統(tǒng)， 剛-柔-熱耦合， 改進的中心差分法

引言

航天器太陽帆板或者太陽帆板的組成構件多為板式結構.在軌航天器沿軌道飛行或者航天器的太陽帆板展開過程中，太陽帆板會受到周期性、不均勻的太陽光照射，太陽帆板上下表面的熱流密度的差異引起太陽帆板沿厚度方向溫度梯度，導致沿厚度方向的熱膨脹不均勻，從而誘發(fā)太陽帆板的熱振動，在一定程度上會影響航天器的姿態(tài)運動.

Thornton和Kim［1］在不考慮航天器運動的情況下建立了卷筒式太陽電池陣熱與結構耦合的動力學模型，研究了熱與結構耦合的動力學特性，然后Johnston和Thornton［2］在考慮航天器的姿態(tài)運動的情況下，研究了在軌運動航天器在太陽輻射熱流作用下太陽帆板的熱誘發(fā)振動對航天器姿態(tài)運動的影響.在此基礎上，Johnston和Thornton［3］進一步開展了動力學實驗，驗證了理論結果的準確性.筆者［4］將熱與結構耦合的動力學模型推廣到薄板結構，研究了剛-柔-熱耦合特性.以上研究工作均基于線彈性理論，僅適用于小變形問題.

由于溫度上升幅值較大，需要考慮太陽帆板熱膨脹引起的軸向變形對彎曲變形的影響.此外，由于太陽帆板具有質量輕和尺度大的特點，彈性變形較為顯著，因此，為了保證計算結果的準確性，在建立動力學方程的過程中需要考慮幾何非線性項.在幾何非線性動力學研究方面，國內外學者針對大變形板殼構件的力學研究開展了一系列工作.Oguamanam等學者［5］針對熱沖擊作用下中心剛體和復合材料層合殼組成的多體系統(tǒng)，利用Von Karman應變與位移的關系建立了考慮幾何非線性的剛-柔耦合動力學方程，在研究剛-柔耦合特性的同時，研究了幾何非線性項對動力學特性的影響.筆者［6］從Green應變與位移的關系式出發(fā)，建立了熱載荷作用下作任意大范圍運動柔性梁的幾何非線性動力學方程，對溫度升高引起的動力軟化問題進行了分析.以上研究工作均假定溫度場是給定的，無法對剛-柔-熱耦合機理進行分析.此外，在推導彎曲應變的過程中利用了板殼曲率的近似表達式，對于大變形問題，容易導致數(shù)值計算誤差.

目前，Shabana提出的絕對節(jié)點坐標法已廣泛應用于大變形柔性多體系統(tǒng)動力學分析，該方法的特點是選取柔性梁或柔性板各單元節(jié)點相對慣性基的位移和斜率作為廣義坐標，建立動力學方程［7-8］，這樣就不需要引入大范圍運動變量和變量位移坐標，廣義坐標全部是全局坐標，得到的質量陣為常值陣，廣義力的表達式也比較簡單.Dombrowski等學者［9］用絕對節(jié)點坐標法對大變形薄板的動力學建模理論進行研究，建立了矩形薄板的二維有限元模型.Mikkola［10］等學者進一步考慮剪切變形，建立了三維矩形薄板的有限元模型.為了驗證絕對節(jié)點坐標法的正確性，Yoo等學者［11-12］首先用快速相機對大變形薄板進行實驗研究，并將測量技術推廣到大變形復雜結構多體系統(tǒng).趙將等學者［13］針對簡化的“IKAROS”自旋展開太陽帆系統(tǒng)，采用結合自然坐標方法與絕對節(jié)點坐標方法的絕對坐標方法對其進行建模，研究了黏彈性太陽帆薄膜自旋展開過程的動力學特性.

為了保證剛-柔耦合動力學仿真精度，本文采用絕對節(jié)點坐標法對熱載荷作用下中心剛體-薄板多體系統(tǒng)進行動力學建模.不計薄板的面外剪切變形，基于Kirchhoff假設，從格林應變和曲率與絕對位移的精確關系式出發(fā)，推導了非線性廣義彈性力陣，用絕對節(jié)點坐標法建立了大變形復合材料薄板的有限元離散的動力學變分方程.在此基礎上進一步考慮剛體姿態(tài)運動、彈性變形和溫度變化的相互耦合作用，建立了熱流密度與絕對節(jié)點坐標之間的關系式，推導了中心剛體-矩形板的剛-柔-熱三者耦合的熱傳導方程和動力學方程.考慮到Runge-kutta法的計算效率較低，本文用廣義-α法求解動力學方程，用改進的中心差分法求解熱傳導方程，差分后的方程通過牛頓迭代法耦合求解.通過仿真算例對無耦合模型，剛-柔-熱耦合的仿真結果進行分析，得到了對工程實際具有參考價值的結論.

1 熱載荷作用下矩形薄板的動力學方程

大變形矩形薄板如圖1所示，設板的長度為a，寬度為b，厚度為h，體密度為ρ.建立固結在地面的慣性坐標系有限元方法將矩形薄板等分為n＝n1×n2個單元，對第e個單元建立單元坐標系則板單元的長度和寬度分別為ae＝a／n1，be＝b／n2.設（x，y，0）為薄板中面上任意一點在該單元坐標系下的坐標，基于絕對節(jié)點坐標法，中面上對應點k0的絕對位置矢量在慣性坐標系下的坐標陣為：

圖1 大變形矩形薄板Fig.1 Thin rectangular plate with large deformation

其中，S（x，y）為形函數(shù)矩陣，由文獻［11］給出.qe（t）為單元節(jié)點的絕對位置坐標，位置坐標分別對x，y的一階導數(shù)，以及位置坐標對x和y的兩階導數(shù)在慣性基下的坐標陣，qe（t）表達式分別為：

其中，qke（t），k＝1，…，4為單元4個節(jié)點的絕對節(jié)點坐標列陣，具體表述方式如下：

設q（t）為總體絕對節(jié)點坐標列陣，qe（t）與q（t）的關系為qe（t）＝Beq（t），其中，Be為單元e的布爾陣.

基于Kirchhoff假設，不計板的剪切變形，為了表示任意一點k的應變，引入點k所在中面上相對應點k0的應變列陣ε0，設n0（x，y，t）為變形后沿法線的單位矢量，基于Kirchhoff假設，法線n0（x，y，t）與中面保持垂直.格林應變的表達式可寫成：

其中，ε0為中面的面內應變列陣，可表示為：

κ為中面的曲率列陣，可表示為：

1.1 質量陣和廣義外力陣的推導

針對薄板問題，可以忽略絕對位移沿單元坐標系Ze方向的變化，板單元e的慣性力的虛功率為：

將式（2）代入上式：

其中，單元質量陣Me為：

矩形薄板的慣性力的虛功率為：

利用關系式（2），板單元體力的虛功率為：

其中，單元體力列陣Qge為：

體積力的虛功率為：

1.2 廣義彈性力陣的推導

對中面的面內應變求變分，并利用關系式（2），得到：

其中，

基于Kirchhoff假設，中面的單位法向量在慣性基下的坐標陣n0可表示為：

其中，n為rx與ry的叉積，n為矢量n的模，可表示為：

對（19）求變分，利用式（2），得到：

將式（18）代入（7），得到κ的表達式為：

對式（24）求變分，并利用關系式（2），得到：

下面推導薄板在熱載荷作用下，當溫度變化ΔT時，板結構的熱應變可表示為［6］：

薄板結構的應力與應變的本構關系可表示為：

其中，E為彈性模量，γ為泊松比.

板單元彈性力的虛功率為：

將式（5）和（28）代入式（31）可得板彈性力的虛功率的表達式為：

化簡后得到板的彈性力虛功表達式：

其中，單元廣義彈性力陣的具體表達式為：

根據(jù)虛功原理，動力學變分方程為：

1.3 中心剛體和薄板多體系統(tǒng)的動力學方程

設中心剛體的質量為mr，關于過質心主軸的轉動慣量陣為J′，以中心剛體質心C為原點建立連體基質心C相對慣性基的位置矢量在慣性基下的坐標陣為rC，連體基相對慣性基的角速度矢量在連體基下的坐標陣為ω′C，中心剛體的動力學變分方程為：

其中，F(xiàn)C為作用于中心剛體的外力的主矢在慣性基下的坐標陣，M′C為作用于中心剛體的外力關于質心的主矩在連體基下的坐標陣.

中心剛體的動力學變分方程為：

中心剛體和矩形薄板多體系統(tǒng)的動力學變分方程為：

設柔性多體系統(tǒng)的約束方程為：

柔性多體系統(tǒng)含拉格朗日乘子的第一類拉格朗日動力學方程為：

其中，

Φq為約束方程的雅可比矩陣，λ為由拉格朗日乘子構成的列陣.

將式（45）和（46）聯(lián)立構成微分代數(shù)混合方程：

2 熱傳導方程

2.1 熱傳導方程

在三維情況下，溫度場變量T在笛卡爾坐標系中滿足的熱傳導變分方程為：

其中，c是材料的比熱容，k是材料的熱傳導系數(shù)，Q是物體內部的熱源密度，q+、q-分別為上下表面的熱流密度分別為上下表面的外界環(huán)境溫度，hc為換熱系數(shù).其余邊界處絕熱，僅有上下表面是傳熱的.

設TN與TM分別是板中面上任意一點的溫度變化和沿厚度方向的溫度梯度，T0為參考溫度，則對于比熱容不大的情況，非中面上任意一點的溫度可表示成T＝T0+TN+zTM，即非中面上任意一點的溫度變化為ΔT＝TN+zTM.將公式（51）和（52）代入熱傳導變分方程（50），得到：

將式（54）代入變分方程（53），則有限元離散柔性板的熱傳導變分方程為

2.2 考慮剛-柔-熱耦合的熱流密度q的推導

假設上下表面的換熱系數(shù)為0，板內部的熱源密度為0，本文考慮的熱流密度q主要來源于宇宙中太陽輻射的熱流.圖2所示為太陽帆板上任意點的中面法線矢量和太陽光輻射矢量.

圖2 有大范圍運動的太陽帆板Fig.2 The solar plate with large overallmotion

其中，n為矢量n的模，由式（20）給出.將式（66）代入式（67），太陽直接輻射角系數(shù)可以表示為：

上式表明，太陽直接輻射系數(shù)為矩陣的絕對節(jié)點坐標q的函數(shù)，說明太陽直接輻射角系數(shù)與太陽帆板的剛體運動和彈性變形三者相互耦合.

為了區(qū)分究竟是單元的上表面還是下表面對著太陽，太陽直接輻射角系數(shù)Γα可以修改為：

求得太陽直接輻射角系數(shù)Γ后，考慮剛-柔-熱耦合的熱流密度可表示為：

將式（52）和（54）代入（35），單元廣義彈性力陣Qfe3和Qfe4的具體表達式為

根據(jù)熱流密度的表達式（68）可知，本文建立的熱流密度函數(shù)與薄板的絕對節(jié)點坐標有關系，而薄板的絕對節(jié)點坐標q和剛體運動以及柔性體的彈性變形相關，那么基于公式（48）和（63）建立的動力學模型即為本文提出剛-柔-熱三者耦合模型.若在建模中忽略中心剛體的姿態(tài)運動和薄板的彈性變形對熱流密度影響，則任意點的法線矢量陣n為常值，此時計算得到熱流密度與薄板的絕對節(jié)點坐標q并無關系，僅與太陽常數(shù)S有關，那么基于此種工況的熱流密度建立的動力學模型就是無耦合的動力學模型.由此可見，剛-柔-熱三者耦合的動力學模型是最精確的力學模型.

3 數(shù)值計算方法

3.1 剛-柔-熱耦合多體系統(tǒng)的動力學方程

將熱傳導方程和剛體與板的動力學方程聯(lián)合，結合公式（48）和（65）建立剛-柔-熱耦合多體系統(tǒng)動力學方程如下：

3.2 數(shù)值計算方法

為了在保證計算精度的同時保證剛-柔-熱耦合動力學方程的計算效率，本文采用對動力學方程用廣義-α法構造差分格式，而熱傳導方程用改進的中心差分法構造差分格式，兩者聯(lián)立后用Newton-Raphson方法求解微分代數(shù)混合方程（71）.計算步驟為：

（1）動力學方程運用廣義-α法構造的差分格式如下：

對熱傳導方程采用改進的中心差分格式后可變?yōu)椋?/p>

（2）將差分格式（76）～（78）代入方程（73）后使用牛頓-拉斐遜迭代法求解聯(lián)合方程組如下：

4 數(shù)值仿真

本算例中，用中心剛體模擬衛(wèi)星的艙體，薄板結構模擬太陽帆板.中心剛體關于質心連體基的轉動慣量為Jx＝Jy＝Jz＝200kg·m2，太陽帆板長為7.5m，寬為0.5m，厚度為0.0103m，采用航空航天工程中常用的低密度、低導熱系數(shù)的各向同性材料，其彈性模量E＝1.93×109N／m2，密度36.8kg／m3，導熱系數(shù)k＝1.5W／（m·K），比熱容c＝921J／（kg ·K），熱膨脹系數(shù)α＝2.3×10-5K-1，帆板的表面吸收率αm＝0.79，太陽常數(shù)S＝1350W／m2.太陽帆板與中心剛體通過固支鉸連接.設結構阻尼系數(shù)為c＝0.01.

設太陽入射角為α0，即與慣性基的夾角為π-α0，則在慣性基下的坐標陣為s＝下面對衛(wèi)星太陽帆板姿態(tài)機動的整個驅動過程進行模擬，以研究太陽熱輻射對整個動力學過程的影響.設中心剛體質心C相對慣性基靜止，繞軸作定軸轉動.初始時刻，中心剛體連體基相對慣性基的三個姿態(tài)角均為0，參考溫度為0.本算例在中心剛體上施加方向驅動力偶，其時變歷程如圖3所示.太陽輻射的入射角為α0＝65°.

圖3 驅動力偶Fig.3 Driving torque

圖4 中心剛體C方向轉角Fig.4 Rotational angle of the hub alongCaxis

圖4和圖5分別為計算得到的中心剛體轉角和角速度的時變圖，第一階段中心剛體繞軸加速轉動，當時間t大于25s時，由于驅動力偶為0，停止加速，角速度在平均值附近振動，第二階段角速度的平均值為常數(shù).當時間t大于175s時，開始第三階段的減速.當時間t大于200s，驅動力偶為0，角速度在0附近振動，停在固定角度，也就是第四階段.從圖中可以看出，非耦合模型和剛-柔-熱耦合模型的得到的中心剛體方向轉角基本吻合.

圖5 中心剛體C方向角速度Fig.5 Angular velocity of the hub alongCaxis

圖6為帆板中面上角點P的沿厚度溫度梯度曲線.圖7為帆板中面上角點P沿方向縱向變形.圖8為沿方向橫向變形.從圖中可以看到，對于剛-柔-熱三者耦合模型，由于熱流密度的表達式中考慮了角系數(shù)與薄板絕對位移坐標的關系，當航天器受到驅動力偶作用作定軸轉動時，沿厚度方向的溫度梯度呈現(xiàn)周期性變化，當姿態(tài)機動到位時，沿厚度方向的溫度梯度保持不變.而無耦合模型由于熱流密度為常數(shù)，沿厚度方向的溫度梯度始終保持不變.此外，對于無耦合模型，由于沿厚度方向溫度梯度為常數(shù)，P點橫向變形和沿方向縱向變形均不受剛體轉動的影響.對于剛-柔-熱三者耦合模型，P點橫向變形為高頻和低頻振蕩曲線的迭加，而低頻振蕩的頻率正是沿厚度方向溫度梯度振動的頻率.除了橫向變形，沿方向縱向變形在緩慢增大的同時（熱膨脹引起），也呈現(xiàn)周期性的波動.可以看出，當姿態(tài)機動到位時，本文非線性模型得到的橫向變形平均值的絕對值明顯大于不考慮縱向變形影響的傳統(tǒng)的線性模型.

圖6 P點厚度方向溫度梯度Fig.6 Temperature gradient of PointPin thickness direction

圖7 P點沿C方向縱向變形Fig.7 Longitudinal deformation of PointPalongCaxis

圖8 P點沿C方向橫向變形Fig.8 Transverse deformation of PointPalongCaxis

5 結論

基于熱傳導方程和有限元離散化方法，建立了熱傳導方程的變分方程，考慮剛體姿態(tài)運動、彈性變形和溫度變化的相互耦合，建立了熱流密度與剛體姿態(tài)坐標和彈性坐標的精確關系式.結合剛-柔耦合動力學方程和熱傳導方程.考慮熱邊界條件隨剛體運動和彈性變形的變化，建立了具有通用形式的剛-柔-熱耦合多體系統(tǒng)的動力學方程.為了有效地提高計算效率，將改進的中心差分法和廣義-α法相結合，求解熱傳導方程和動力學方程，差分后的方程通過牛頓迭代法進行求解.

對無耦合模型和剛-柔-熱三者耦合模型的仿真結果進行比較表明，中心剛體的姿態(tài)運動對沿厚度方向溫度梯度和熱變形的影響不容忽視，表現(xiàn)在沿厚度方向溫度梯度呈現(xiàn)周期性的振蕩.此外，橫向變形和縱向變形曲線是低頻振動和高頻振動曲線的迭加，低頻振動頻率等于沿厚度方向溫度梯度振動的頻率.進一步研究還發(fā)現(xiàn)，本文非線性模型得到的橫向變形平均值的絕對值明顯大于不考慮軸向變形對橫向變形影響的傳統(tǒng)的線性模型.從而說明當軸向變形隨著溫度升高逐漸增大時，軸向變形在方向的分量對橫向變形具有一定影響，不容忽視.

1 Thornton E A，Kim Y A.Thermally induced bending vibrations of a flexible rolled-up solar array.Journal of Spacecraft and Rockets，1993，30（4）：438～448

2 Johnston J D，Thornton E A.Thermally induced attitude dynamics ofa spacecraftwith a flexible appendage.Journal of Guidance，Control and Dynamics，1998，21（4）：581～587

3 Johnston JD，Thornton E A.Thermally induced dynamics of satellite solar panels.Journal of Spacecraft and Rockets，2000，37（5）：604～613

4 潘科琪.曲梁和板殼結構多體系統(tǒng)剛-柔耦合動力學研究［博士學位論文］.上海：上海交通大學，2012（Pan K Q.Rigid-flexible coupling dynamics for curved beam and shell structuremultibody systems.［PhD Thesis］Shanghai：Shanghai Jiao Tong University，2012（in Chinese））

5 Oguamanam D CD，Hansen JS，Heppler G R.Nonlinear transient response of thermally loaded laminated panels.ASME Journal of Applied Mechanics，2004，71（1）：49～56

6 Liu JY，Lu H.Thermal effecton the deformation of a flexible beam with large kinematical driven overall motions.European Journal of Mechanics A／Solids，2007，26（1）：137～151

7 Berzeri M，Shabana A A.Development of simple models for the elastic forces in absolute nodal co-ordinate formulation.Journal of Sound and Vibration，2000，235（4）：539～565

8 Berzeri M，Shabana A A.Study of the centrifugal stiffening effect using the finite element absolute nodal coordinate formulation.Multibody System Dynamics，2002，7（4）：357～387

9 Dombrowski SV.Analysis of large flexible body deformation in multibody systems using absolute coordinates.Multibody System Dynamics，2002，8（4）：409～432

10 Mikkola A M，Shabana A A.A non-incremental finite element procedure for the analysis of large deformation of plates and shells inmechanical system applications.Multibody System Dynamics，2003，9（3）：283～309

11 Yoo W S，Kim M S，Mun SH，Sohn JH.Large displacement of beam with basemotion：flexiblemultibody simulations and experiments.Computing Methods Applied Mechanical Engineering，2006，195（50）：7036～7051

12 YooW S，Kim M S，Mun SH，Sohn JH.Developments of multibody system dynamics：computer simulations and experiments.Multibody System Dynamics，2007，18（1）：35～5

13 趙將，劉鋮，田強，胡海巖.黏彈性薄膜太陽帆自旋展開動力學分析.力學學報，2013，45（2）：746～754（Zhao J，Liu C，Tian Q，Hu H Y.Dynamic analysis of spinning deployment of a solar sail composed of viscoelastic membranes.Chinese Journal of Theoretical and Applied Mechanics，2013，45（2）：746～754（in Chinese） ）

RIGID-FLEXIBLE-THERMALCOUPLING DYNAM IC FORMULATION FOR HUB-PLATE MULTIBODY SYSTEM*

Zhang Weihua Liu Jinyang?
（Department of Engineering Mechanics，Shanghai Jiaotong，University，Shanghai200240，China）

In this paper，rigid-flexible-thermal coupling dynamic formulation for hub rectangular platemultibody system with large deformation under thermal load is investigated.Based on Kirchhoff assumption，and the nonlinear relationships among the Green strain，the curvature and the absolute coordinates，the nonlinear elastic force is derived.Additionally，the finite element discretized dynamic variational equations of the rectangular plate are derived using absolute nodal coordinate formulation.In order to investigate the rigid-flexible-thermal coupling effect，the relationship between the heat flux and the absolute coordinates is also obtained.And through the constrained equations of themultibody system，the rigid-flexible-thermal coupling equations are derived，which include the heat conduction equations and the Lagrange dynamic equations of the first kind.Moreover，in order to improve the simulation efficiency，modified central differencemethod and generalized-a method are combined to solve the heat conduction equations and the dynamic equations.And the discretized equations are solved by using Newton-Raphsonmethod.Satellite and platemultibody system applied with solar heat flux is taken as examples to verify the effectiveness of the formulation.Comparison of the results obtained by non-couplingmethod and rigidflexible-thermal couplingmethod shows that the effect of the rigid bodymotion on the temperature gradient in the thickness direction and thermally induced deformation is significant，which should be taken into account.Furthermore，due to the inclusion of the geometric nonlinear terms，the influence of the axial deformation caused by the temperature increase on the transverse deformation is also taken into consideration.

hub-platemultibody system， rigid-flexible-thermal coupling， modified central differencemethod

10.6052／1672-6553-2015-88

2015-4-8收到第1稿，2015-11-12收到修改稿.

*國家自然科學基金資助項目（11272203，11132007）

?通訊作者E-mail：liujy＠sjtu.edu.cn

Received 8 April2015，revised 12 November 2015.

*The project supported by the National Natural Science Foundation of China（11272203，11132007）

?Corresponding author E-mail：liujy＠sjtu.edu.cn