李建宏宇，龐賀偉

（1.北京衛(wèi)星環(huán)境工程研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100094；2.中國(guó)空間技術(shù)研究院，北京 100094）

微振動(dòng)是影響航天器有效載荷特性的重要因素[1]。 大型航天器微振動(dòng)涉及到結(jié)構(gòu)系統(tǒng)、控制系統(tǒng)、光學(xué)系統(tǒng)等。傳統(tǒng)的集中在單系統(tǒng)的設(shè)計(jì)和評(píng)估方法已經(jīng)不能滿足實(shí)際工程的需求，建立結(jié)構(gòu)-控制-光學(xué)一體化模型是目前比較有效的分析手段[2]。21 世紀(jì)初，國(guó)外在SIM[3]、JWST[4]、TPF[5]等高分辨率空間望遠(yuǎn)鏡的研制中開(kāi)發(fā)了IMOS[6]、DOCS[7-8]和IME[9]等一體化建模軟件。國(guó)內(nèi)一些學(xué)者對(duì)一體化建模方法進(jìn)行了介紹和綜述性的研究[10-11]，并應(yīng)用到基于整星級(jí)的微振動(dòng)仿真[12-13]。

目前，結(jié)構(gòu)-控制-光學(xué)一體化建模方法主要應(yīng)用于概念設(shè)計(jì)階段，航天器結(jié)構(gòu)系統(tǒng)相對(duì)簡(jiǎn)單，建立的一體化模型的狀態(tài)空間維數(shù)較低，不存在數(shù)值求解困難的問(wèn)題。如SIM 中一體化模型狀態(tài)空間的維數(shù)為316[3]；JWST 中狀態(tài)空間的維數(shù)為320[4]；國(guó)內(nèi)的整星級(jí)一體化模型狀態(tài)空間的維數(shù)達(dá)到 1×103量級(jí)[12-13]。在實(shí)際大型航天器詳細(xì)設(shè)計(jì)階段，由于結(jié)構(gòu)復(fù)雜，且自由度數(shù)高，一體化模型的狀態(tài)空間維數(shù)將達(dá)到1×104量級(jí)，其計(jì)算成本將增長(zhǎng)到千倍以上。在建立一體化模型計(jì)算微振動(dòng)響應(yīng)時(shí)，需要考慮模型和數(shù)值算法的計(jì)算效率。

文中使用狀態(tài)空間法和附加剛度法建立一體化模型，根據(jù)兩種模型的特點(diǎn)分別選取合適的數(shù)值計(jì)算方法，并制定了微振動(dòng)頻率響應(yīng)計(jì)算策略，極大地提高了微振動(dòng)響應(yīng)的計(jì)算效率。

1 大型航天器微振動(dòng)的分系統(tǒng)模型

1.1 結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)模型

結(jié)構(gòu)動(dòng)力學(xué)模型是基于有限元模型建立的。結(jié)構(gòu)動(dòng)力學(xué)方程為：

由于航天器結(jié)構(gòu)復(fù)雜，自由度數(shù)高，一般利用模態(tài)疊加法對(duì)結(jié)構(gòu)進(jìn)行頻率響應(yīng)分析?？紤]前r 階模態(tài)，其中r 由需要考察的頻率范圍決定。一般情況下，考察前n Hz 的頻率響應(yīng)需要至少提取模態(tài)頻率小于2n Hz 的所有模態(tài)。

有限元模型由航天器結(jié)構(gòu)設(shè)計(jì)單位提供，內(nèi)容包括特征值信息和特定節(jié)點(diǎn)的特征向量信息。一般情況下，不會(huì)要求輸出全部自由度上的特征向量信息。這是由于大型航天器結(jié)構(gòu)有限元模型的自由度達(dá)到百萬(wàn)以上，導(dǎo)出全部自由度的特征向量信息存在著存儲(chǔ)成本過(guò)大和數(shù)據(jù)存儲(chǔ)和讀取時(shí)間過(guò)長(zhǎng)的問(wèn)題。

以模態(tài)坐標(biāo)及其一階導(dǎo)數(shù)作為線性系統(tǒng)的狀態(tài)量，則結(jié)構(gòu)動(dòng)力學(xué)的線性控制系統(tǒng)為：

線性控制系統(tǒng)的輸入分為兩部分：擾動(dòng)力（力矩）和控制系統(tǒng)的控制力（力矩）。線性控制系統(tǒng)的輸出分為兩部分，分別是控制系統(tǒng)中傳感器自由度上的位移（轉(zhuǎn)動(dòng)）和光學(xué)系統(tǒng)中光學(xué)結(jié)構(gòu)自由度上的位移（平動(dòng)和轉(zhuǎn)動(dòng)）。

1.2 控制系統(tǒng)模型

控制系統(tǒng)模型由姿態(tài)控制設(shè)計(jì)單位以線性控制系統(tǒng)的狀態(tài)空間形式提供。該系統(tǒng)的輸入為控制系統(tǒng)中傳感器自由度上的位移（轉(zhuǎn)動(dòng)），輸出為控制系統(tǒng)的控制力（力矩）。

其中輸入 xs與結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)模型中的狀態(tài)量的關(guān)系為：

控制系統(tǒng)在頻域的傳遞函數(shù)為：

1.3 光學(xué)系統(tǒng)模型

光學(xué)系統(tǒng)模型由光學(xué)設(shè)計(jì)單位提供。在微振動(dòng)的力學(xué)環(huán)境下，光學(xué)鏡面之間的相對(duì)運(yùn)動(dòng)為微小運(yùn)動(dòng)，可以用光學(xué)系統(tǒng)的一階差分計(jì)算出光學(xué)靈敏度矩陣。

式中： yo是光軸繞三個(gè)坐標(biāo)軸的轉(zhuǎn)角；為光學(xué)靈敏度矩陣； xo為光學(xué)系統(tǒng)局部模型中光學(xué)結(jié)構(gòu)自由度上的位移（平動(dòng)和轉(zhuǎn)動(dòng)）。

1.4 選擇矩陣

控制模型和光學(xué)模型都是建立在局部模型中。與結(jié)構(gòu)動(dòng)力學(xué)模型建立一體化模型時(shí)，需要將局部的自由度擴(kuò)展到整個(gè)結(jié)構(gòu)的全局自由度上。

設(shè)一個(gè)自由度集合 P = { p1, p2, … , ps}。自由度 pj在局部模型中為第個(gè)自由度，在全局模型中為第個(gè)自由度。局部坐標(biāo)系的總自由度數(shù)為d，全局坐標(biāo)系的總自由度數(shù)為n。則選擇矩陣為n ×d 矩陣，其中第行第列的元素為1，其他元素均為0?？梢酝ㄟ^(guò)選擇矩陣 βP將自由度集合P 在局部模型和全局模型之間相互轉(zhuǎn)換。

在一體化模型中，選擇矩陣與特征值矩陣是以βTΦ 的形式成對(duì)出現(xiàn)。根據(jù)選擇矩陣的定義可得出，對(duì)于自由度集合P：

2 大型航天器結(jié)構(gòu)-控制-光學(xué)一體化模型

2.1 通過(guò)狀態(tài)空間法建立一體化模型

將結(jié)構(gòu)動(dòng)力學(xué)模型和控制系統(tǒng)模型整合，建立新的狀態(tài)方程：

令一體化模型的狀態(tài)量為動(dòng)力學(xué)系統(tǒng)中的模態(tài)坐標(biāo)及其一階導(dǎo)數(shù)和控制系統(tǒng)中的狀態(tài)量：

則有：

將光學(xué)系統(tǒng)的輸出作為一體化模型的輸出：

則一體化模型在頻域的傳遞函數(shù)為：

光軸三個(gè)方向的轉(zhuǎn)動(dòng)歐拉角與擾動(dòng)力（力矩）在頻域的關(guān)系為：

這是通過(guò)狀態(tài)空間法建立一體化模型的微振動(dòng)響應(yīng)計(jì)算公式。

尾礦的組成較為復(fù)雜，我國(guó)尾礦主要成分有石英、長(zhǎng)石、石灰石、滑石、白云石、云母、高嶺石、石榴石和綠泥石等。按照尾礦的主要組成成分，尾礦可分為以下5類(lèi)：

2.2 通過(guò)附加剛度建立一體化模型

在控制系統(tǒng)模型中，通過(guò)傳遞函數(shù)計(jì)算控制系統(tǒng)的控制力（力矩）與控制系統(tǒng)中傳感器自由度上的位移（轉(zhuǎn)動(dòng)）的關(guān)系為：

通過(guò)選擇矩陣將局部模型中控制系統(tǒng)的控制力（力矩）與傳感器自由度上的位移（轉(zhuǎn)動(dòng)）擴(kuò)展到全局模型中。在全局模型中，控制力（力矩）與位移（轉(zhuǎn)動(dòng)）的關(guān)系為：

將控制力（力矩）作為附加剛度加入到動(dòng)力學(xué)方程中，動(dòng)力學(xué)方程變?yōu)椋?/p>

使用模態(tài)疊加法，將動(dòng)力學(xué)方程縮聚到Φ 構(gòu)成的r 階線性子空間后，轉(zhuǎn)換到頻域：

對(duì)于光學(xué)系統(tǒng)模型，光軸三個(gè)方向的轉(zhuǎn)動(dòng)歐拉角與模態(tài)坐標(biāo)的關(guān)系為：

這是通過(guò)附加剛度法建立一體化模型的微振動(dòng)響應(yīng)計(jì)算公式。

2.3 微振動(dòng)響應(yīng)求解過(guò)程中的數(shù)值方法

對(duì)于兩種建模方法計(jì)算微振動(dòng)響應(yīng)，主要的計(jì)算成本均為求解線性方程組。其基本形式為：

對(duì)于狀態(tài)空間法建立的一體化模型， A = ( jω E-；對(duì)于附加剛度建立的一體化模型，。在兩種方法中，矩陣A 不屬于對(duì)稱(chēng)矩陣、三角矩陣、埃爾米特（Hermitian）矩陣、黑森貝格（Hessenberg）矩陣中的任意一種矩陣，可采用LU 矩陣分解法進(jìn)行求解。

對(duì)于狀態(tài)空間法建立的一體化模型，A 為稀疏矩陣，可以使用稀疏矩陣的LU 矩陣分解法進(jìn)行進(jìn)一步加速；對(duì)于附加剛度建立的一體化模型，A 由兩部分組成，其中一部分為對(duì)角矩陣，可以采用迭代法進(jìn)行求解。

以一個(gè)典型的大型航天器一體化模型為例。狀態(tài)空間法建立的一體化模型中，方陣A 的維數(shù)為16 702，根據(jù)估算[14]，使用LU 分解的計(jì)算復(fù)雜度在1×1012量級(jí)；方陣A 中非零元素個(gè)數(shù)為1×106量級(jí)，根據(jù)估算[15]，使用稀疏矩陣LU 分解的計(jì)算復(fù)雜度在1×109量級(jí)。

附加剛度建立的一體化模型中，方陣A 的維數(shù)為8348，根據(jù)估算，使用LU 分解的計(jì)算復(fù)雜度在1×1011量級(jí)。使用迭代法進(jìn)行求解線性方程組的步驟為：令；選取 x0=D 作為初始近似 值；迭 代 公 式 為；當(dāng)?shù)螖?shù)超過(guò)設(shè)置的最大迭代次數(shù)或者時(shí)，迭代結(jié)束。

迭代過(guò)程中，使用矩陣乘法的結(jié)合律，可以減少計(jì)算復(fù)雜度。迭代過(guò)程中的計(jì)算為矩陣相乘。經(jīng)過(guò)估算，迭代一次的計(jì)算復(fù)雜度在1×106量級(jí)。迭代初始近似值的物理意義是無(wú)控制系統(tǒng)下系統(tǒng)的頻率響應(yīng)。由于控制系統(tǒng)僅在低頻起到抑制振動(dòng)的效果，而對(duì)高頻振動(dòng)影響不大。從物理角度分析，在低頻段需要迭代次數(shù)較多，而在高頻段需要迭代次數(shù)較少。

迭代收斂的必要條件是矩陣G 的任意一種范數(shù)小于1。在整個(gè)頻段的計(jì)算中，當(dāng)頻率逼近0 或者頻率逼近固有頻率時(shí)，可能存在某些頻率點(diǎn)使用該迭代公式無(wú)法收斂的情況。

選取了5 個(gè)典型的頻率點(diǎn)對(duì)兩種建模方法使用四種數(shù)值計(jì)算方法進(jìn)行微振動(dòng)響應(yīng)計(jì)算。在典型的大型航天器微振動(dòng)預(yù)示計(jì)算中，需要計(jì)算的頻率點(diǎn)在1×103量級(jí)。以2800 頻率點(diǎn)為例，預(yù)估四種數(shù)值計(jì)算方法在全頻段的時(shí)間成本，結(jié)果見(jiàn)表1。

表1 典型頻率點(diǎn)計(jì)算時(shí)間成本Tab.1 Time costs on typical frequency point calculation

通過(guò)表1 可以驗(yàn)證估算的計(jì)算復(fù)雜度的正確性。使用迭代法求解附加剛度建立的一體化模型微振動(dòng)響應(yīng)的計(jì)算效率遠(yuǎn)高于其他方法，但在某些頻率點(diǎn)可能出現(xiàn)迭代發(fā)散的情況。使用稀疏矩陣LU 分解求解狀態(tài)空間建立的一體化模型微振動(dòng)響應(yīng)的計(jì)算效率雖然不是最高，但仍然在可接受范圍，且適用于所有頻率點(diǎn)。

結(jié)合兩種算法的特點(diǎn)，可以制定全頻段微振動(dòng)響應(yīng)計(jì)算策略：首先使用迭代法求解附加剛度建立的一體化模型微振動(dòng)響應(yīng)，當(dāng)在某些頻率點(diǎn)，迭代法發(fā)散或者收斂困難時(shí)，使用稀疏矩陣LU 分解求解狀態(tài)空間建立的一體化模型微振動(dòng)響應(yīng)。這樣既可保證整個(gè)頻段微振動(dòng)響應(yīng)的計(jì)算效率，又可避免特殊頻率點(diǎn)求解失敗。

3 數(shù)值算例

3.1 某大型航天器一體化建模與微振動(dòng)預(yù)示

某大型航天器有限元模型由400 000+個(gè)節(jié)點(diǎn)和400 000+個(gè)單元組成，自由度數(shù)為2 000 000+。利用Nastran 軟件對(duì)有限元模型進(jìn)行模態(tài)分析，取前8348階（小于400 Hz 的所有模態(tài)）特征值，并輸出特征向量在擾動(dòng)力、控制力、傳感器、光學(xué)系統(tǒng)節(jié)點(diǎn)自由度上的分量。

航天器的主要擾源有CMG、制冷機(jī)等。擾動(dòng)輸入為150 個(gè)自由度上的擾動(dòng)力和擾動(dòng)力矩。擾動(dòng)力和擾動(dòng)力矩由諧波疊加組成。控制模型的控制輸入為3個(gè)轉(zhuǎn)動(dòng)自由度，控制輸出為3 個(gè)控制力矩。光學(xué)模型的輸入為5 個(gè)節(jié)點(diǎn)共計(jì)30 個(gè)平動(dòng)自由度和轉(zhuǎn)動(dòng)自由度，輸出為光軸三個(gè)方向的轉(zhuǎn)動(dòng)歐拉角。

分別通過(guò)狀態(tài)空間法和附加剛度法建立一體化模型?？疾祛l率從0.02～200 Hz 的微振動(dòng)頻率響應(yīng)。在0.02～20 Hz 的低頻段，取1000 個(gè)頻率點(diǎn)（間隔0.02 Hz）；在20～200 Hz 的中高頻段取1800 個(gè)頻率點(diǎn)（間隔0.1 Hz）。根據(jù)2.3 小節(jié)中的計(jì)算策略求解整個(gè)頻段共計(jì)2800 個(gè)頻率點(diǎn)的微振動(dòng)頻率響應(yīng)。

從表1 可知，使用迭代法求解附加剛度法的一體化模型中，迭代法一次的計(jì)算成本是稀疏矩陣LU 分解法求解狀態(tài)空間法計(jì)算成本的1/1000 左右。據(jù)此設(shè)置最大迭代次數(shù)設(shè)為100，設(shè)置迭代收斂的收斂精度 ε = 10-4x0，保證迭代法求解附加剛度法的一體化模型與稀疏矩陣LU 分解法求解狀態(tài)空間法的相對(duì)誤差在0.01%左右。

3.2 計(jì)算結(jié)果

低頻段和中高頻段的光軸繞X 軸轉(zhuǎn)動(dòng)角的頻率響應(yīng)如圖1 所示，同時(shí)繪制了無(wú)控制系統(tǒng)下的微振動(dòng)頻率響應(yīng)。在低頻段，相比于無(wú)控制系統(tǒng)，有控制系統(tǒng)明顯抑制了微振動(dòng)響應(yīng)；在中高頻段，控制系統(tǒng)的影響減弱，有控制系統(tǒng)的頻率響應(yīng)和無(wú)控制的頻率響應(yīng)幾乎重合。

圖1 光軸繞X 軸轉(zhuǎn)動(dòng)角的頻率響應(yīng)Fig.1 Frequency response of rotation angle of optical axis around X axis

在整個(gè)頻段共2800 個(gè)頻率點(diǎn)中，共有6 個(gè)頻率點(diǎn)使用迭代法求解附加剛度建立的一體化模型微振動(dòng)響應(yīng)存在困難，其中4 個(gè)頻率點(diǎn)（0.02、0.04、0.06、0.10 Hz）接近0 Hz，2 個(gè)頻率點(diǎn)（0.58 Hz 接近第8階模態(tài)0.5807 Hz，1.12 Hz 接近第10 階頻率1.1163 Hz）逼近固有頻率?？偤臅r(shí)10.35 s，其中使用迭代法計(jì)算2794 個(gè)頻率點(diǎn)耗時(shí)7.20 s，使用稀疏LU 分解法計(jì)算6 個(gè)頻率點(diǎn)耗時(shí)3.15 s。微振動(dòng)的計(jì)算成本得到了很好的控制，使未來(lái)針對(duì)微振動(dòng)響應(yīng)的優(yōu)化提供了保障。

迭代法的迭代次數(shù)如圖2 所示?？梢钥闯?，在小于40 Hz 以下的中低頻段，迭代次數(shù)較高；在大于40 Hz 以上的中高頻段，迭代次數(shù)較低。這說(shuō)明了迭代初始值在中高頻段與最終迭代值相近，其物理意義為無(wú)控制系統(tǒng)的頻率響應(yīng)與有控制系統(tǒng)的頻率響應(yīng)幾乎相等。即從數(shù)學(xué)角度和物理角度相互印證了控制系統(tǒng)在大型航天器微振動(dòng)抑制中對(duì)中高頻段的效果較弱的結(jié)論。

圖2 數(shù)值計(jì)算迭代次數(shù)Fig.2 Number of iterations for numerical computation

4 結(jié)語(yǔ)

大型航天器由于結(jié)構(gòu)復(fù)雜自由度數(shù)高，在建立結(jié)構(gòu)-控制-光學(xué)一體化模型求解微振動(dòng)響應(yīng)時(shí)，需要考慮計(jì)算成本。文中使用狀態(tài)空間法和附加剛度法建立一體化模型并推導(dǎo)出微振動(dòng)響應(yīng)計(jì)算公式，并根據(jù)計(jì)算公式的特性分別采用稀疏矩陣LU 分解法和迭代法兩種數(shù)值計(jì)算方法。使用迭代法求解附加剛度法的一體化模型效率最高，但在特殊頻率點(diǎn)存在迭代發(fā)散的現(xiàn)象，這時(shí)可以換用稀疏矩陣LU 分解法求解狀態(tài)空間法的一體化模型。數(shù)值算例證明了方法的高效性。此方法是適用于實(shí)際工程中大型航天器微振動(dòng)預(yù)示的高效算法。