李志華，李廣，沈漢武，樊志華

杭州電子科技大學(xué) 機(jī)械工程學(xué)院，杭州 310018

具有剛性、非線性和強(qiáng)耦合特點(diǎn)的動(dòng)力學(xué)系統(tǒng)廣泛存在于機(jī)械、土木、航空航天等領(lǐng)域。例如航天器系統(tǒng)，由于存在柔性件振動(dòng)變形以及大范圍運(yùn)動(dòng)位移和小范圍振動(dòng)位移之間的相互耦合，其動(dòng)力學(xué)模型就呈現(xiàn)剛性、非線性和強(qiáng)耦合的特點(diǎn)[1-2]。對(duì)這類剛性系統(tǒng)進(jìn)行仿真求解十分困難，而對(duì)其動(dòng)力學(xué)模型的精確高效求解是對(duì)其進(jìn)行精確控制的基礎(chǔ)[3]。目前數(shù)值求解的方法主要包括直接數(shù)值求解法和降階求解法。

直接數(shù)值求解法以Wilson-θ和Newmark為主[4-6]，該方法需要調(diào)用迭代算法，計(jì)算成本高、仿真效率低，且仿真精度的提高存在局限性。降階求解法通過將動(dòng)力學(xué)方程轉(zhuǎn)化為低階的常微分方程，然后采用傳統(tǒng)的數(shù)值積分方法進(jìn)行求解，相比于直接數(shù)值求解法，它在仿真效率、計(jì)算成本以及仿真精度上均存在一定的優(yōu)勢(shì)。

降階求解法所使用的傳統(tǒng)數(shù)值積分方法有Euler法、Adams多步法、龍格庫塔法和Gear法等。Euler和龍格庫塔法屬于單步法，Gear和Adams法屬于多步法。其中，Euler法是最簡(jiǎn)單的數(shù)值方法，它在求解剛性問題時(shí)，需要顯著地減少積分步長(zhǎng)，以保持穩(wěn)定性，但這會(huì)導(dǎo)致誤差的積累，不利于剛性問題的求解；Gear法初始值較難確定，求解高振蕩方程時(shí)會(huì)失效，同時(shí)方程的階數(shù)一般較大[7]；隱式的龍格庫塔方法求解剛性問題時(shí)穩(wěn)定性好、精度高，但是s階的龍格庫塔法在求解一個(gè)r維的剛性問題時(shí)，每步都需要求解一個(gè)由s×r個(gè)方程聯(lián)立的方程組，在實(shí)際工程應(yīng)用中，其計(jì)算量過于龐大，導(dǎo)致實(shí)際求解時(shí)的效率不高[8]。相比于龍格庫塔方法，Adams多步法的優(yōu)點(diǎn)是計(jì)算格式簡(jiǎn)單、每步計(jì)算量相對(duì)較小，同時(shí)也具備較高的求解精度，缺點(diǎn)是使用隱式Adams多步法在求解剛性問題時(shí)，相比于顯式算法仍會(huì)有較大的計(jì)算量。由此可見，單純地使用傳統(tǒng)的隱式算法或顯式算法來求解動(dòng)力學(xué)中的剛性問題，都存在著一定的困難。

量化狀態(tài)系統(tǒng)(Quantized State System，QSS)算法是一種新的數(shù)值積分方法，由Zeigler和Lee[9]首次提出，并由 Kofman和Junco[10]首次算法實(shí)現(xiàn)。與基于時(shí)間離散的傳統(tǒng)方法不同，QSS將所有的狀態(tài)變量離散化，以“量子”代替“步長(zhǎng)”，通過計(jì)算所有狀態(tài)變量每次變化一個(gè)量子所需要的最短時(shí)間，來進(jìn)行下一步積分的推進(jìn)。在求解非剛性問題時(shí)，QSS算法穩(wěn)定性強(qiáng)，同時(shí)計(jì)算過程不需要迭代，因此仿真效率得到了極大的提高。QSS的狀態(tài)變量軌跡是分段線性的，其求解精度相對(duì)于傳統(tǒng)算法并不占優(yōu)勢(shì)。為了提高其算法精度，QSS2和QSS3采用高次曲線來對(duì)狀態(tài)變量變化的時(shí)間進(jìn)行估計(jì)，利用狀態(tài)變量的高階導(dǎo)數(shù)來改進(jìn)近似值[11-13]。在求解非剛性問題時(shí)，QSS2和QSS3的求解精度得到了提高，要好于QSS。但在求解剛性問題時(shí)，QSS、QSS2、QSS3均存在可能的數(shù)值振蕩現(xiàn)象；此外，在求解剛性問題時(shí)，QSS2和QSS3的求解效率要遠(yuǎn)低于QSS[12-14]。

為了解決QSS求解剛性問題時(shí)所出現(xiàn)的數(shù)值振蕩現(xiàn)象，Migoni和Kofman[15]提出了BQSS(Backward QSS)算法，但該算法在求解非線性剛性問題時(shí)，由于誤差范圍較大，可能會(huì)出現(xiàn)“假性平衡”現(xiàn)象，影響仿真的推進(jìn)。后續(xù)Migoni提出了一階LIQSS(Linearly Implicit QSS)算法，它將量化變量q看作狀態(tài)變量x的預(yù)測(cè)值，避免了某些擾動(dòng)項(xiàng)對(duì)求解剛性問題所造成的影響。但是LIQSS算法要求剛性問題中的雅可比矩陣的對(duì)角線要存在較大的元素，否則可能會(huì)出現(xiàn)“假性振蕩”現(xiàn)象，會(huì)對(duì)仿真求解的精度和穩(wěn)定性造成較大的影響[16]。之后，Migoni等[17]又提出了改進(jìn)的LIQSS算法，但在求解剛性較大的微分方程時(shí)，仍然存在求解效率較低或者無法求解的缺陷。國(guó)內(nèi)學(xué)者目前研究還較少，只有文獻(xiàn)[18-20]對(duì)QSS算法進(jìn)行了初步應(yīng)用。

本文將QSS算法和隱式多步法相結(jié)合，提出一種自適應(yīng)多步校正算法。通過對(duì)柔性航天器動(dòng)力學(xué)的仿真求解，驗(yàn)證了算法的可行性。從仿真精度和效率2方面，將該算法與幾種典型的傳統(tǒng)算法(ODE23tb、ODE15s、ODE45等)以及QSS算法進(jìn)行了性能對(duì)比。

1 量化狀態(tài)系統(tǒng)

常微分方程表示的狀態(tài)方程系統(tǒng)為

(1)

式中:x(t)∈Rn為系統(tǒng)的狀態(tài)向量；u(t)為輸入向量，QSS算法將式(1)方程重新定義為

(2)

式中:q(t)為系統(tǒng)的量化變量，每個(gè)量化變量qi(t)均通過式(3)的遲滯量化函數(shù)得到，

(3)

式中:ΔQi為量子；qi(t-)為上一次計(jì)算的量化變量。

圖1 QSS算法的流程圖

2 自適應(yīng)多步校正算法

2.1 AMCQSS

本文提出的自適應(yīng)多步校正算法(Adaptive Multi-step Correction algorithm based on QSS，AMCQSS)旨在有效提高仿真精度的同時(shí)保證仿真的效率，并且拓展QSS算法求解剛性問題的應(yīng)用范圍。本算法以QSS算法為基礎(chǔ)、同時(shí)結(jié)合了隱式多步法的思想：QSS算法作為顯式算法，具有較好的全局誤差控制特性，且其求解過程不需要迭代，保證了仿真求解的效率；多步法能夠充分運(yùn)用前面已經(jīng)求得的多個(gè)點(diǎn)的信息來校正目前所需求解的值；同時(shí)本算法還可以根據(jù)計(jì)算過程中求得的導(dǎo)數(shù)差值來自行選擇二步法或三步法，以獲得更高的求解精度；此外本算法將量化變量值作為狀態(tài)變量的預(yù)測(cè)值，來扼制剛性系統(tǒng)中可能出現(xiàn)的仿真數(shù)值振蕩，以提高算法的穩(wěn)定性和仿真求解的精度。

AMCQSS方法將式(1)近似為

(4)

式中:q為量化變量，這里作為求x的一階導(dǎo)數(shù)的預(yù)測(cè)值。

(5)

(6)

此時(shí)，狀態(tài)變量仿真推進(jìn)時(shí)間為

(7)

式中：k為仿真執(zhí)行的步數(shù)。

當(dāng)k≤2時(shí)，狀態(tài)變量值的計(jì)算公式為

(8)

當(dāng)k>2時(shí)，AMCQSS算法使用多步法思想對(duì)第k步的狀態(tài)變量進(jìn)行校正，本算法中融入了隱式多步法中的二步法和三步法思想。在求解剛性問題時(shí)，系統(tǒng)方程的曲線可能會(huì)在某一點(diǎn)突變，此時(shí)的曲線斜率會(huì)劇烈變化，該點(diǎn)的導(dǎo)數(shù)與前幾點(diǎn)的導(dǎo)數(shù)相比較會(huì)產(chǎn)生較大的變化，這時(shí)使用二步法或三步法思想來對(duì)該點(diǎn)的狀態(tài)變量進(jìn)行校正，其校正求得的2個(gè)狀態(tài)變量在數(shù)值上可能差距較大。為了保證求解的精度，將二步法與三步法求得的導(dǎo)數(shù)分別和第k步導(dǎo)數(shù)相比，選擇導(dǎo)數(shù)值與第k步導(dǎo)數(shù)值相差較小的多步法。此時(shí)所選擇的多步法校正后得到的第k步狀態(tài)變量值是位于另一種多步法校正后得到的狀態(tài)變量值和校正前的狀態(tài)變量值之間，無論該點(diǎn)真實(shí)情況偏向于哪一邊，均不會(huì)出現(xiàn)偏差過大的情況，這樣可以有效地控制誤差，避免誤差積累影響到之后的仿真推進(jìn)。導(dǎo)數(shù)差計(jì)算公式為

(9)

根據(jù)差值選擇狀態(tài)變量的計(jì)算公式

(10)

(11)

(12)

然后按式(7)～式(10)計(jì)算，最終系統(tǒng)每次推進(jìn)的仿真時(shí)間為

Δt=min(Δtj)

(13)

2.2 AMCQSS算法的實(shí)現(xiàn)

算法實(shí)現(xiàn)(以偽代碼的形式)如下：

While(t

for

until (|xj-qj=ΔQj|)

qj=xj+ΔQj

qj=xj-ΔQj

if (k≤2) then

else

else

end if

end if

else

if (k≤2)then

else

else

end if

end if

end for

Δt=min(Δtj) ∥(系統(tǒng)每推進(jìn)一步，仿真的時(shí)間取狀態(tài)變量的最小躍遷時(shí)間)

end while

2.3 算法的收斂性與穩(wěn)定性分析

2.3.1 收斂性

(14)

假設(shè)λ(t)和λ1(t)分別為式(4)和式(14)的初始解，且λ(0)=λ1(0)，則

f(λ(τ),u(τ))]dτ

(15)

根據(jù)Euclidean范數(shù)，可得

(16)

設(shè)定M為函數(shù)f的Lipschitz常數(shù)，則式(16)有

MtΔx

(17)

其中λ(t),λi(t)均連續(xù)，且M為正數(shù)，則根據(jù)Gronwall-Bellman不等式將式(17)改為

(18)

其中M和t不依賴于Δx，則

(19)

由式(19)可得，AMCQSS算法具有收斂性。

2.3.2 穩(wěn)定性

假設(shè)式(1)是線性時(shí)不變系統(tǒng)，可將式(1)改寫為

(20)

使用AMCQSS算法可將式(20)寫成

(21)

式中：A為Hurwitz矩陣；B為系統(tǒng)矩陣，則AMCQSS的誤差為

|e(t)|≤|V|·|Re(Λ)-1·Λ|·|V-1|·ΔQ

(22)

式中：Λ為矩陣A的特征值的對(duì)角矩陣；V為矩陣A的特征向量矩陣；Re為取復(fù)數(shù)的實(shí)數(shù)部分。

由式(22)可得，AMCQSS算法的全局誤差總是有界，因此AMCQSS算法具有穩(wěn)定性。

3 柔性航天器的動(dòng)力學(xué)模型

圖2是柔性航天器中的柔性多體衛(wèi)星力學(xué)模型，將衛(wèi)星本體(星體)假設(shè)為無變形的中心剛體、忽略鉸連接處的摩擦力、將太陽帆板視為固連在星體上的一塊均質(zhì)薄板且相對(duì)星體無轉(zhuǎn)動(dòng)[23-26]。

圖2 柔性多體衛(wèi)星系統(tǒng)力學(xué)模型

圖2中，B為星體，Ai表示為與星體連接的第i個(gè)撓性附件。引入與軌道參考坐標(biāo)系等同的慣性坐標(biāo)系Oxyz、星體固連坐標(biāo)系Oxbybzb和撓性附件固連坐標(biāo)系Pxiyizi，忽略衛(wèi)星自身運(yùn)動(dòng)相對(duì)標(biāo)稱位置的小擾動(dòng)[25]。Ob為星體質(zhì)心；Pi為撓性附件與星體的連接點(diǎn)；Rb,k為點(diǎn)O到星體質(zhì)量元dmb,k的矢徑；rb,k為點(diǎn)Ob到質(zhì)量元dmb,k的矢徑；Rai,j為點(diǎn)O到質(zhì)量元dmai,j的矢徑；xO為星體質(zhì)心相對(duì)于點(diǎn)O的小位移攝動(dòng)；rpi為質(zhì)心點(diǎn)Ob和點(diǎn)Pi之間的矢徑；rai,j為Pi到撓性附件質(zhì)量元dmai,j上的矢徑；δai,j為dmai,j的振動(dòng)變形矢量[25-27]。

使用模態(tài)綜合—混合坐標(biāo)法以及Euler-Lagrange方程對(duì)柔性多體衛(wèi)星進(jìn)行動(dòng)力學(xué)建模，并將其轉(zhuǎn)化成常微分方程形式(由于篇幅所限，本文不進(jìn)行詳細(xì)推導(dǎo)，具體過程可參考文獻(xiàn)[25-26])：

(23)

式中:ΛL和ΛR分別為左右帆板的模態(tài)剛度陣；FsaiL和FsaiR分別為整星轉(zhuǎn)動(dòng)時(shí)左右帆板的剛?cè)狁詈详?；FtaiL和FtaiR分別為整星平動(dòng)時(shí)左右帆板的剛?cè)狁詈详?。為了方便?jì)算，記FsaiL1為FsaiL矩陣的第1行，F(xiàn)saiL2為FsaiL矩陣的第2行，以此類推。M為衛(wèi)星質(zhì)量矩陣；As為衛(wèi)星的姿態(tài)變換陣(歐拉角轉(zhuǎn)動(dòng)順序?yàn)?-1-3)：

其中：cos表示為c；sin表示為s；α為繞慣性系y軸旋轉(zhuǎn)α角，得到過渡坐標(biāo)系x1y1z1；θ為繞過渡坐標(biāo)系x1y1z1中x1軸旋轉(zhuǎn)θ角，得到過渡坐標(biāo)系x2y2z2；φ為繞過渡坐標(biāo)系x2y2z2中z2軸旋轉(zhuǎn)φ角，得到星體固連坐標(biāo)系xbybzb，完成慣性系到星體固連坐標(biāo)系的轉(zhuǎn)換。

4 仿 真

4.1 仿真背景

ΛR=ΛL=diag{900，14 400，16 900，250 000，280 900，360 000}，左右帆板的剛?cè)狁詈详?單位為kg1/2·m)為

FtaiR=

FtaiL=

⑤ 系統(tǒng)各個(gè)狀態(tài)變量的相對(duì)誤差計(jì)算公式為

(24)

式中:um[k]為各個(gè)算法求得的狀態(tài)變量值；udassl[k]為DASSL求解器在誤差設(shè)定為10-9的情況下求得的狀態(tài)變量值，在此作為基準(zhǔn)值[28]。

4.2 仿真結(jié)果對(duì)比分析

對(duì)柔性多體衛(wèi)星動(dòng)力學(xué)模型進(jìn)行仿真求解，得到圖3～圖5所示的衛(wèi)星角度及角速度變化圖，其中QSS和ODE45算法無解。從圖3～圖5中可以看出：ODE23t出現(xiàn)了明顯的發(fā)散，說明ODE23t不適用于該剛性問題的求解；ODE23tb雖然最終結(jié)果趨于穩(wěn)定，但中間過程上下波動(dòng)較大，其求解剛性問題的穩(wěn)定性較差；AMCQSS和ODE15s 2種算法波動(dòng)均很小，且整個(gè)過程收斂性和穩(wěn)定性均較好。由表1可以看出，其中求解結(jié)果最好的傳統(tǒng)算法是ODE15s；相比于ODE15s，AMCQSS算法在求解效率上提高了1倍多，在求解精度上提高了4倍多。

圖3 衛(wèi)星α角轉(zhuǎn)角變化圖

圖4 衛(wèi)星θ角轉(zhuǎn)角變化圖

圖5 衛(wèi)星φ角轉(zhuǎn)角速度

表1 柔性多體衛(wèi)星動(dòng)力學(xué)方程的仿真結(jié)果

5 結(jié) 論

針對(duì)具有剛性、非線性和強(qiáng)耦合特點(diǎn)的動(dòng)力學(xué)求解難題，本文提出了一種有別于傳統(tǒng)方法的自適應(yīng)多步校正算法AMCQSS，通過對(duì)柔性航天器動(dòng)力學(xué)的仿真求解，得到以下結(jié)論：

1) ODE45和QSS是顯式算法，在求解剛性問題時(shí)，會(huì)出現(xiàn)無解情況，很難滿足實(shí)際工程中剛性問題的求解需要。

2) ODE23t和ODE23tb算法穩(wěn)定性較差，不適于強(qiáng)剛性問題的求解。

3) ODE15s和AMCQSS算法穩(wěn)定性和收斂性均較好，但AMCQSS算法的求解精度和效率更高，算法性能優(yōu)于QSS算法和ODE15s等傳統(tǒng)算法。