本發(fā)明屬于空間結(jié)構(gòu)動力學(xué)計算技術(shù)領(lǐng)域，具體涉及一種空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)動力學(xué)建模與計算方法。

背景技術(shù)：

近年來，我國航天科技發(fā)展迫切需要掌握柔性空間結(jié)構(gòu)在軌展開技術(shù)，以滿足衛(wèi)星通信、天基對地觀測和深空探測等重大需求。這類空間結(jié)構(gòu)一般主要由剛性桿、柔性梁、板、殼和柔性繩索等構(gòu)件組成，各構(gòu)件間除鉸節(jié)點處的位置約束外，還存在大量同步、滑移等非線性約束。同時這類空間結(jié)構(gòu)展開尺度大、柔軟部件多、構(gòu)形復(fù)雜，其展開過程會呈現(xiàn)系統(tǒng)大范圍運動與柔軟部件大變形之間的非線性耦合動力學(xué)。美國、俄羅斯等國在其航天任務(wù)中，曾發(fā)生多起結(jié)構(gòu)展開失敗。由于這類結(jié)構(gòu)的展開動力學(xué)地面實驗難度大、無法完全抵消重力影響，故其展開過程的動力學(xué)數(shù)值模擬是確保結(jié)構(gòu)展開成功的關(guān)鍵技術(shù)。

針對此類空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)，經(jīng)典多體系統(tǒng)動力學(xué)方法的研究思路是:在柔性結(jié)構(gòu)體上附著浮動坐標(biāo)系,以其相對于絕對坐標(biāo)系的運動來描述柔體大范圍運動,通過模態(tài)縮聚近似描述柔性結(jié)構(gòu)體相對浮動坐標(biāo)系的變形。該類方法一般引入小轉(zhuǎn)動、小變形假設(shè),將柔性結(jié)構(gòu)體的彈性變形與剛體位移分開考慮，進行解耦處理。通常這種方法只能靠大量增加單元數(shù)量來處理柔性結(jié)構(gòu)體的大變形、大位移問題,從而導(dǎo)致計算規(guī)模大、計算效率低。此外,這種方法不能精確地反映空間復(fù)雜柔性結(jié)構(gòu)展開過程中的力學(xué)本質(zhì)特性。由于空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)存在很多非線性約束、運動副間隙及復(fù)雜的載荷條件,目前商業(yè)軟件在系統(tǒng)動力學(xué)建模、數(shù)值求解效率等方面存在許多困難。

技術(shù)實現(xiàn)要素：

有鑒于此，本發(fā)明提供了一種空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)動力學(xué)建模與計算方法，為柔性結(jié)構(gòu)多體系統(tǒng)提供了一種簡單有效的動力學(xué)建模方法，提高了計算效率。

實現(xiàn)本發(fā)明的技術(shù)方案如下：

步驟101：建三維模型和賦予參數(shù)：根據(jù)所要求解的實物，建立該實物的三維模型，輸入所建模型的材料參數(shù)、幾何參數(shù)以及空間結(jié)構(gòu)拓撲關(guān)系；

步驟102：網(wǎng)格劃分：基于步驟101賦予參數(shù)后的三維模型進行網(wǎng)格劃分，利用絕對節(jié)點坐標(biāo)法、幾何精確法或等幾何分析方法確定網(wǎng)格單元類型，將劃分后的三維模型定義為網(wǎng)格模型；

步驟103：根據(jù)步驟102劃分的網(wǎng)格確定各個網(wǎng)格單元的初始節(jié)點位置和初始節(jié)點速度，基于單元初始節(jié)點位置和初始節(jié)點速度計算網(wǎng)格單元廣義坐標(biāo)向量q和廣義速度向量

步驟104：根據(jù)步驟103獲得的廣義坐標(biāo)向量和廣義速度向量，設(shè)置網(wǎng)格模型的外部載荷及邊界條件，基于網(wǎng)格模型的邊界條件，建立網(wǎng)格模型的約束方程，并設(shè)置仿真總時長t以及仿真時間步長h；

步驟105：根據(jù)步驟102-104獲得的網(wǎng)格單元類型、廣義坐標(biāo)向量、廣義速度向量和約束方程，基于第一類拉格朗日方程建立多體系統(tǒng)動力學(xué)方程，如下：

其中m為系統(tǒng)質(zhì)量陣，為廣義加速度向量，f(·)為系統(tǒng)的彈性力向量，q(·)為系統(tǒng)的廣義外力向量，φ(·)和φq(·)分別為步驟104獲得的約束方程向量和約束方程向量對廣義坐標(biāo)向量的偏導(dǎo)數(shù)矩陣，(·)q為函數(shù)對廣義坐標(biāo)向量的偏導(dǎo)數(shù)，λ為拉格朗日乘子向量，t為仿真時間,t∈[0,t]，t為矩陣轉(zhuǎn)置；

步驟106：對所建的多體系統(tǒng)動力學(xué)方程進行時域離散，獲得時域離散后的多體系統(tǒng)動力學(xué)方程，對離散后的多體系統(tǒng)動力學(xué)方程進行迭代，并對迭代過程中生成的線性方程組進行并行求解，輸出廣義坐標(biāo)向量q以及拉格朗日乘子向量λ，至此完成對空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)動力學(xué)問題的建模與求解。

進一步地，所述步驟106具體為：

步驟601：根據(jù)步驟105得到多體系統(tǒng)動力學(xué)方程，對多體系統(tǒng)動力學(xué)方程進行離散，離散后多體系統(tǒng)動力學(xué)方程見公式(2)；

步驟602：對離散后的多體系統(tǒng)動力學(xué)方程進行求解得到第n+1步的廣義坐標(biāo)向量qn+1、廣義速度向量和拉格朗日乘子向量λn+1；

步驟603：根據(jù)步驟602判斷時間tn+1是否小于總仿真時間t，若當(dāng)前時間tn+1不小于總仿真時間t，則計算停止；否則，n自動加1，重復(fù)執(zhí)行步驟601-603，計算廣義坐標(biāo)向量qn+1、廣義速度向量令時間tn+1＝tn+h，直到tn+1大于或等于t，輸出廣義坐標(biāo)向量qn+1以及拉格朗日乘子向量λn+1，其中qn+1為q，λn+1為λ。

進一步地，對多體系統(tǒng)動力學(xué)方程進行離散采用如下方法：

應(yīng)用廣義-alpha隱式時間積分算法對多體系統(tǒng)動力學(xué)方程進行如下離散：

其中

式(4)中，β和γ為決定計算精度與效率的算法參數(shù)，取γ≥1/2和β≥(1/2+γ)²/4。n為迭代次數(shù)，h為仿真時間步長，為了耗散系統(tǒng)高頻響應(yīng)，引入新的算法矢量參數(shù)a，其中矢量參數(shù)a滿足關(guān)系：

式(4)和式(5)中各參數(shù)的選取方法如下:

其中為算法的譜半徑，αm和αf均為廣義-alpha隱式時間積分算法參數(shù)。

進一步地，本步驟中用q均代表qn+1，λ均代表λn+1，所述步驟602進一步包括：

步驟6021：設(shè)置求解精度tol；

步驟6022：采用newton-rapson迭代求解所述公式(2)，獲得如下線性代數(shù)方程組：

其中δq和δλ是當(dāng)前時刻t系統(tǒng)廣義坐標(biāo)向量和拉格朗日乘子向量的增量，其中

和為滿足如下關(guān)系的算法參數(shù):

其中，和均為中間變量，ψ和γ均為中間矩陣，為函數(shù)對廣義速度向量的偏導(dǎo)數(shù)；

步驟6023：根據(jù)步驟6021獲得的線性方程組計算當(dāng)前殘差向量的模r；

其中，||·||為模；

步驟6024：根據(jù)步驟6022對線性方程組的求解更新當(dāng)前廣義坐標(biāo)向量以及拉格朗日乘子向量；

步驟6025：根據(jù)步驟6023獲得的殘差向量模判斷迭代結(jié)果是否收斂，若r>tol，則根據(jù)步驟6024更新獲得的廣義坐標(biāo)向量以及拉格朗日乘子向量，重復(fù)執(zhí)行步驟6022-6025，直至r<tol或r＝tol，輸出步驟6024獲得的廣義坐標(biāo)向量q以及拉格朗日乘子向量λ。

有益效果：

1)本發(fā)明采用絕對節(jié)點坐標(biāo)法、幾何精確梁理論和等幾何分析等非線性有限元法對空間復(fù)雜柔性結(jié)構(gòu)進行離散，充分考慮了柔性結(jié)構(gòu)的運動與變形耦合問題，為空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)的動力學(xué)建模提供了手段。

2)本發(fā)明采用區(qū)域分解和并行計算技術(shù)，大幅度提高了空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)動力學(xué)仿真計算效率。

附圖說明

圖1為對含復(fù)雜多柔體系統(tǒng)動力學(xué)問題進行建模與計算方法流程圖。

圖2為某桁架索網(wǎng)天線結(jié)構(gòu)示意圖。

圖3為本發(fā)明所開發(fā)的部分可選的單元類型示意圖。

圖4為軸承軸瓦模型結(jié)構(gòu)示意圖。

圖5為同步齒輪模型結(jié)構(gòu)示意圖。

圖6為移動副模型結(jié)構(gòu)示意圖。

圖7為對系統(tǒng)動力學(xué)方程高效數(shù)值積分求解算法流程圖。

圖8為對系統(tǒng)離散非線性方程高效求解算法流程圖。

圖9為本發(fā)明建模和求解算法方法應(yīng)用于某桁架索網(wǎng)天線，并進行數(shù)值仿真天線展開部分時刻結(jié)構(gòu)示意圖。

具體實施方式

下面結(jié)合附圖并舉實施例，對本發(fā)明進行詳細描述。

對于空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)，通常由眾多柔性桿件和柔性索網(wǎng)組成，在展開過程中這些構(gòu)件會產(chǎn)生很大變形和大范圍轉(zhuǎn)動，此時基于小變形和小轉(zhuǎn)動假設(shè)的多體建模分析方法已不能滿足工程需求。對于此類柔性構(gòu)件的建模主要采用非線性有限元法：如絕對節(jié)點坐標(biāo)法、幾何精確法和等幾何分析方法等。其中shabana等基于連續(xù)介質(zhì)力學(xué)和有限元法提出了絕對節(jié)點坐標(biāo)法(ancf)。與經(jīng)典多體系統(tǒng)動力學(xué)方法相比，該方法在慣性坐標(biāo)系下采用統(tǒng)一的插值函數(shù)描述柔體的剛體運動與彈性變形，非常適用于描述同時具有大轉(zhuǎn)動和大變形的柔體動力學(xué)。此外，幾何精確梁理論是simo等以有限應(yīng)變梁理論為基礎(chǔ)發(fā)展的一種非線性有限元方法，該法能夠高效并準(zhǔn)確地處理大變形與大轉(zhuǎn)動問題，近年來也被廣泛應(yīng)用于空間可展開結(jié)構(gòu)的動力學(xué)仿真。為了統(tǒng)一幾何模型與力學(xué)分析模型，實現(xiàn)兩者的無縫集成，在有限元領(lǐng)域，hughes等提出的一種新的以樣條理論為基礎(chǔ)的數(shù)值計算方法：等幾何分析(iga)。

綜上所述，在空間復(fù)雜柔性結(jié)構(gòu)多體系統(tǒng)動力學(xué)展開分析階段，如何考慮并解決柔性構(gòu)件變形與大范圍運動的耦合問題，并提高數(shù)值仿真計算效率成為航天科技發(fā)展亟待解決的技術(shù)問題。本發(fā)明提出了一種對含大變形構(gòu)件的空間結(jié)構(gòu)多體系統(tǒng)動力學(xué)建模與計算方法。

如圖1所示，實施典型案例為對如圖2所示的某環(huán)形桁架索網(wǎng)天線機構(gòu)展開動力學(xué)建模與計算方法。

該方法包括如下步驟：

步驟101：建三維模型和賦予參數(shù)：根據(jù)所要求解的實物，建立該實物的三維模型，輸入所建模型的材料參數(shù)、幾何參數(shù)以及空間結(jié)構(gòu)拓撲關(guān)系。具體包括：根據(jù)環(huán)形桁架索網(wǎng)天線的幾何參數(shù)以及空間結(jié)構(gòu)拓撲關(guān)系，在proe軟件中建立環(huán)形桁架索網(wǎng)天線三維模型，并輸入環(huán)形桁架天線不同部件的材料參數(shù)。由于金屬絲反射網(wǎng)密度很小，且附著在前張力網(wǎng)上，因而忽略金屬絲反射網(wǎng)的影響。

步驟102：網(wǎng)格劃分：基于步驟101賦予參數(shù)后的三維模型進行網(wǎng)格劃分，利用絕對節(jié)點坐標(biāo)法、幾何精確法或等幾何分析方法確定網(wǎng)格單元類型，將劃分后的三維模型定義為網(wǎng)格模型；根據(jù)模型性質(zhì)以及精度要求，選取合適的單元類型，如圖3所示，剛體模型有：基于自然坐標(biāo)方法描述的剛體單元、基于自然坐標(biāo)方法描述的縮減剛體單元，基于絕對節(jié)點坐標(biāo)方法描述的參考節(jié)點剛體單元；柔體單元有：基于絕對節(jié)點坐標(biāo)方法描述的全參數(shù)梁單元、基于絕對節(jié)點坐標(biāo)方法描述的縮減梁單元、基于絕對節(jié)點坐標(biāo)方法描述的斜率不連續(xù)的梁單元、基于絕對節(jié)點坐標(biāo)描述的矩形薄板單元以及基于絕對節(jié)點坐標(biāo)描述的三角形薄板單元等，其中絕對節(jié)點坐標(biāo)法、幾何精確法或等幾何分析方法均具有解決耦合特性的特點。

步驟103：確定模型廣義坐標(biāo)向量和廣義速度向量：根據(jù)步驟102劃分的網(wǎng)格確定各個網(wǎng)格單元的初始節(jié)點位置和初始節(jié)點速度，基于單元初始節(jié)點位置和初始節(jié)點速度計算網(wǎng)格單元廣義坐標(biāo)向量q和廣義速度向量通過將環(huán)形桁架索網(wǎng)天線收攏到一定角度，從而確定索網(wǎng)和桁架初始廣義坐標(biāo)和廣義速度。

步驟104：根據(jù)步驟103獲得的廣義坐標(biāo)向量和廣義速度向量，設(shè)置網(wǎng)格模型的外部載荷及邊界條件，基于網(wǎng)格模型的邊界條件，建立網(wǎng)格模型的約束方程，并設(shè)置仿真總時長t以及仿真時間步長h。該環(huán)形桁架索網(wǎng)天線仿真總時長為500s，仿真時間步長設(shè)為1e-4s。其中外部載荷為索網(wǎng)的預(yù)張力，約束為索網(wǎng)和桁架之間的鉸接、桁架桿件之間的運動副以及桁架與伸展臂之間的固接。下面簡要介紹三種主要運動副的建模：旋轉(zhuǎn)鉸約束、同步齒輪約束和移動副約束。

其中關(guān)節(jié)旋轉(zhuǎn)鉸約束建立方式如下：對旋轉(zhuǎn)的關(guān)節(jié)處建立剛性的軸承軸瓦模型，如圖4所示。剛性的軸承與軸瓦與相連接的構(gòu)件之間為固定連接，則相連構(gòu)件的旋轉(zhuǎn)可以通過軸承軸瓦的相對轉(zhuǎn)動來描述。采用自然坐標(biāo)描述的軸承軸瓦之間的旋轉(zhuǎn)鉸約束方程可表示為：

其中φ¹為旋轉(zhuǎn)鉸約束，分別由位移約束φ1和旋轉(zhuǎn)角約束φ2組成。和分別為全局坐標(biāo)系下軸承和軸瓦端點位置坐標(biāo)。e¹和e²分別是局部坐標(biāo)系下軸瓦和軸承的單位向量，θ為軸承繞軸瓦的轉(zhuǎn)動角度，其中

其中l(wèi)b和lj分別是軸承和軸瓦的長度，u^b和v^b是位于軸瓦端點b的單位向量，u^j和v^j是位于軸承端點j的單位向量。

同步齒輪約束建立方式如下：

如圖5所示，同步齒輪約束保證相鄰兩個面的展開角度一致，即第一跨的橫桿(ab桿)與豎桿之間的夾角和第二跨的橫桿(cd桿)與豎桿之間的夾角相等，其具體的約束方程見下式:

其中φ²為同步齒輪約束，和分別是豎桿與兩橫桿的軸線方向矢量。

如圖6所示，移動副約束建立如下：移動副約束需保證細桿gh能在粗桿ef中自由滑動，同時細桿與粗桿在桁架展開任意位置時均在四邊形abcd內(nèi)。其具體的約束方程見下式：

其中φ³為移動副約束。φ3保證e點與a點重合，c點與g點重合，其中r^a、r^e、r^c和r^g分別為點a、e、c和g在全局坐標(biāo)系下的位置坐標(biāo)。φ4保證兩斜桿能相互滑動，其中和分別為點e、h、f和g沿斜桿方向矢量，l為斜桿長度。φ4中的最后一項是為了保證兩斜桿在平面abcd內(nèi)，這樣式φ3與φ4就構(gòu)成了移動副的約束方程，至此，φ¹、φ²和φ³構(gòu)成網(wǎng)絡(luò)模型的約束方程向量。

步驟105：根據(jù)步驟102-104獲得的網(wǎng)格單元類型、廣義坐標(biāo)向量、廣義速度向量和約束方程，基于第一類拉格朗日方程建立多體系統(tǒng)動力學(xué)方程，如下：

其中m為系統(tǒng)質(zhì)量陣，為廣義加速度向量，是廣義坐標(biāo)向量對時間的二階導(dǎo)數(shù)，f(·)為系統(tǒng)的彈性力向量，q(·)為系統(tǒng)的廣義外力向量，φ(·)和φq(·)分別為步驟104獲得的網(wǎng)絡(luò)模型的約束方程向量和約束方程向量對廣義坐標(biāo)向量的偏導(dǎo)數(shù)矩陣，(·)q為函數(shù)對廣義坐標(biāo)向量的偏導(dǎo)數(shù)，λ為拉格朗日乘子向量，t為仿真時間,t∈[0,t]，t為矩陣轉(zhuǎn)置。

步驟106：利用廣義-alpha隱式時間積分算法，對所建的多體系統(tǒng)動力學(xué)方程進行時域離散，獲得時域離散后的多體系統(tǒng)動力學(xué)方程，采用newton-rapson迭代對離散后的多體系統(tǒng)動力學(xué)方程進行迭代，并對迭代過程中生成的線性方程組進行并行求解，輸出廣義坐標(biāo)向量q以及拉格朗日乘子向量λ。

如圖7所示，步驟601-603描述了所述步驟106，具體為：

步驟601：根據(jù)步驟105得到多體系統(tǒng)動力學(xué)方程，應(yīng)用廣義-alpha隱式時間積分算法進行如下離散：

其中

式(7)中，β和γ為決定計算精度與效率的算法參數(shù)，取γ≥1/2和β≥(1/2+γ)²/4。n為迭代次數(shù)，h為時間積分步長，為了耗散系統(tǒng)高頻響應(yīng)，引入新的算法矢量參數(shù)a，其中矢量參數(shù)a滿足關(guān)系：

式(7)和式(8)中各參數(shù)的選取方法如下:

其中為算法的譜半徑，決定著算法能量耗散分布的頻率范圍，通常取αm和αf均為廣義-alpha隱式時間積分算法參數(shù)。

步驟602：將式(7)帶入式(6)得到離散形式的非線性方程組，并對非線性方程組進行求解得到廣義坐標(biāo)向量qn+1、廣義速度向量和拉格朗日乘子向量λ；

步驟603：根據(jù)步驟602判斷時間tn+1是否小于總仿真時間t，若當(dāng)前時間tn+1不小于總仿真時間t，則計算停止；否則，n自動加1，令時間tn+1＝tn+h重復(fù)執(zhí)行步驟601-603，計算廣義坐標(biāo)向量qn+1、廣義速度向量直到tn+1大于或等于t。輸出廣義坐標(biāo)向量qn+1以及拉格朗日乘子向量λn+1，其中qn+1為q，λn+1為λ。

如圖8所示，步驟6021-6025進一步地描述了步驟602，具體為：

步驟6021：設(shè)置求解精度tol；

步驟6022：采用newton-rapson迭代求解公式(6)，獲得如下線性代數(shù)方程組：

其中δq和δλ是當(dāng)前時刻t系統(tǒng)廣義坐標(biāo)向量和拉格朗日乘子向量的增量，其中

和為滿足如下關(guān)系的算法參數(shù):

其中，和均為中間變量，ψ和γ均為中間矩陣，為函數(shù)對廣義速度向量的偏導(dǎo)數(shù)；

步驟6023：根據(jù)步驟6021獲得的線性方程組計算當(dāng)前殘差向量的模r；

其中，||·||為模；

步驟6024：根據(jù)步驟6022對線性方程組的求解更新當(dāng)前廣義坐標(biāo)向量以及拉格朗日乘子向量；

步驟6025：根據(jù)步驟6023獲得的殘差向量模判斷迭代結(jié)果是否收斂，若r>tol，則根據(jù)步驟6024更新獲得的廣義坐標(biāo)向量以及拉格朗日乘子向量，重復(fù)執(zhí)行步驟6022-6025，直至r<tol或r＝tol，輸出步驟6024獲得的廣義坐標(biāo)向量q以及拉格朗日乘子向量λ，并存儲。至此完成對某環(huán)形桁架索網(wǎng)天線機構(gòu)展開動力學(xué)建模與高效計算。圖9給出了某環(huán)形桁架天線展開過程中6個指定時刻的系統(tǒng)構(gòu)型，該圖明顯展示了環(huán)形桁架天線展開過程中的不同步現(xiàn)象相應(yīng)的地面試驗結(jié)果驗證了該數(shù)值仿真結(jié)果的有效性。

綜上所述，以上僅為本發(fā)明的較佳實施例而已，并非用于限定本發(fā)明的保護范圍。凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進等，均應(yīng)包含在本發(fā)明的保護范圍之內(nèi)。