大型航天器結構-控制-光學一體化建模與微振動響應快速計算方法

李建宏宇，龐賀偉

（1.北京衛星環境工程研究所 可靠性與環境工程技術重點實驗室，北京 100094；2.中國空間技術研究院，北京 100094）

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

目前，結構-控制-光學一體化建模方法主要應用于概念設計階段，航天器結構系統相對簡單，建立的一體化模型的狀態空間維數較低，不存在數值求解困難的問題。如SIM 中一體化模型狀態空間的維數為316[3]；JWST 中狀態空間的維數為320[4]；國內的整星級一體化模型狀態空間的維數達到 1×103量級[12-13]。在實際大型航天器詳細設計階段，由于結構復雜，且自由度數高，一體化模型的狀態空間維數將達到1×104量級，其計算成本將增長到千倍以上。在建立一體化模型計算微振動響應時，需要考慮模型和數值算法的計算效率。

文中使用狀態空間法和附加剛度法建立一體化模型，根據兩種模型的特點分別選取合適的數值計算方法，并制定了微振動頻率響應計算策略，極大地提高了微振動響應的計算效率。

1 大型航天器微振動的分系統模型

1.1 結構動力學系統模型

結構動力學模型是基于有限元模型建立的。結構動力學方程為：

由于航天器結構復雜，自由度數高，一般利用模態疊加法對結構進行頻率響應分析。考慮前r 階模態，其中r 由需要考察的頻率范圍決定。一般情況下，考察前n Hz 的頻率響應需要至少提取模態頻率小于2n Hz 的所有模態。

有限元模型由航天器結構設計單位提供，內容包括特征值信息和特定節點的特征向量信息。一般情況下，不會要求輸出全部自由度上的特征向量信息。這是由于大型航天器結構有限元模型的自由度達到百萬以上，導出全部自由度的特征向量信息存在著存儲成本過大和數據存儲和讀取時間過長的問題。

以模態坐標及其一階導數作為線性系統的狀態量，則結構動力學的線性控制系統為：

線性控制系統的輸入分為兩部分：擾動力（力矩）和控制系統的控制力（力矩）。線性控制系統的輸出分為兩部分，分別是控制系統中傳感器自由度上的位移（轉動）和光學系統中光學結構自由度上的位移（平動和轉動）。

1.2 控制系統模型

控制系統模型由姿態控制設計單位以線性控制系統的狀態空間形式提供。該系統的輸入為控制系統中傳感器自由度上的位移（轉動），輸出為控制系統的控制力（力矩）。

其中輸入 xs與結構動力學系統模型中的狀態量的關系為：

控制系統在頻域的傳遞函數為：

1.3 光學系統模型

光學系統模型由光學設計單位提供。在微振動的力學環境下，光學鏡面之間的相對運動為微小運動，可以用光學系統的一階差分計算出光學靈敏度矩陣。

式中： yo是光軸繞三個坐標軸的轉角；為光學靈敏度矩陣； xo為光學系統局部模型中光學結構自由度上的位移（平動和轉動）。

1.4 選擇矩陣

控制模型和光學模型都是建立在局部模型中。與結構動力學模型建立一體化模型時，需要將局部的自由度擴展到整個結構的全局自由度上。

設一個自由度集合 P = { p1, p2, … , ps}。自由度 pj在局部模型中為第個自由度，在全局模型中為第個自由度。局部坐標系的總自由度數為d，全局坐標系的總自由度數為n。則選擇矩陣為n ×d 矩陣，其中第行第列的元素為1，其他元素均為0。可以通過選擇矩陣 βP將自由度集合P 在局部模型和全局模型之間相互轉換。

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

2 大型航天器結構-控制-光學一體化模型

2.1 通過狀態空間法建立一體化模型

將結構動力學模型和控制系統模型整合，建立新的狀態方程：

令一體化模型的狀態量為動力學系統中的模態坐標及其一階導數和控制系統中的狀態量：

則有：

將光學系統的輸出作為一體化模型的輸出：

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

光軸三個方向的轉動歐拉角與擾動力（力矩）在頻域的關系為：

這是通過狀態空間法建立一體化模型的微振動響應計算公式。

尾礦的組成較為復雜，我國尾礦主要成分有石英、長石、石灰石、滑石、白云石、云母、高嶺石、石榴石和綠泥石等。按照尾礦的主要組成成分，尾礦可分為以下5類：

2.2 通過附加剛度建立一體化模型

在控制系統模型中，通過傳遞函數計算控制系統的控制力（力矩）與控制系統中傳感器自由度上的位移（轉動）的關系為：

通過選擇矩陣將局部模型中控制系統的控制力（力矩）與傳感器自由度上的位移（轉動）擴展到全局模型中。在全局模型中，控制力（力矩）與位移（轉動）的關系為：

將控制力（力矩）作為附加剛度加入到動力學方程中，動力學方程變為：

使用模態疊加法，將動力學方程縮聚到Φ 構成的r 階線性子空間后，轉換到頻域：

對于光學系統模型，光軸三個方向的轉動歐拉角與模態坐標的關系為：

這是通過附加剛度法建立一體化模型的微振動響應計算公式。

2.3 微振動響應求解過程中的數值方法

對于兩種建模方法計算微振動響應，主要的計算成本均為求解線性方程組。其基本形式為：

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

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

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

附加剛度建立的一體化模型中，方陣A 的維數為8348，根據估算，使用LU 分解的計算復雜度在1×1011量級。使用迭代法進行求解線性方程組的步驟為：令；選取 x0=D 作為初始近似 值；迭 代 公 式 為；當迭代次數超過設置的最大迭代次數或者時，迭代結束。

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

迭代收斂的必要條件是矩陣G 的任意一種范數小于1。在整個頻段的計算中，當頻率逼近0 或者頻率逼近固有頻率時，可能存在某些頻率點使用該迭代公式無法收斂的情況。

選取了5 個典型的頻率點對兩種建模方法使用四種數值計算方法進行微振動響應計算。在典型的大型航天器微振動預示計算中，需要計算的頻率點在1×103量級。以2800 頻率點為例，預估四種數值計算方法在全頻段的時間成本，結果見表1。

表1 典型頻率點計算時間成本Tab.1 Time costs on typical frequency point calculation

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

結合兩種算法的特點，可以制定全頻段微振動響應計算策略：首先使用迭代法求解附加剛度建立的一體化模型微振動響應，當在某些頻率點，迭代法發散或者收斂困難時，使用稀疏矩陣LU 分解求解狀態空間建立的一體化模型微振動響應。這樣既可保證整個頻段微振動響應的計算效率，又可避免特殊頻率點求解失敗。

3 數值算例

3.1 某大型航天器一體化建模與微振動預示

某大型航天器有限元模型由400 000+個節點和400 000+個單元組成，自由度數為2 000 000+。利用Nastran 軟件對有限元模型進行模態分析，取前8348階（小于400 Hz 的所有模態）特征值，并輸出特征向量在擾動力、控制力、傳感器、光學系統節點自由度上的分量。

航天器的主要擾源有CMG、制冷機等。擾動輸入為150 個自由度上的擾動力和擾動力矩。擾動力和擾動力矩由諧波疊加組成。控制模型的控制輸入為3個轉動自由度，控制輸出為3 個控制力矩。光學模型的輸入為5 個節點共計30 個平動自由度和轉動自由度，輸出為光軸三個方向的轉動歐拉角。

分別通過狀態空間法和附加剛度法建立一體化模型。考察頻率從0.02～200 Hz 的微振動頻率響應。在0.02～20 Hz 的低頻段，取1000 個頻率點（間隔0.02 Hz）；在20～200 Hz 的中高頻段取1800 個頻率點（間隔0.1 Hz）。根據2.3 小節中的計算策略求解整個頻段共計2800 個頻率點的微振動頻率響應。

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

3.2 計算結果

低頻段和中高頻段的光軸繞X 軸轉動角的頻率響應如圖1 所示，同時繪制了無控制系統下的微振動頻率響應。在低頻段，相比于無控制系統，有控制系統明顯抑制了微振動響應；在中高頻段，控制系統的影響減弱，有控制系統的頻率響應和無控制的頻率響應幾乎重合。

圖1 光軸繞X 軸轉動角的頻率響應Fig.1 Frequency response of rotation angle of optical axis around X axis

在整個頻段共2800 個頻率點中，共有6 個頻率點使用迭代法求解附加剛度建立的一體化模型微振動響應存在困難，其中4 個頻率點（0.02、0.04、0.06、0.10 Hz）接近0 Hz，2 個頻率點（0.58 Hz 接近第8階模態0.5807 Hz，1.12 Hz 接近第10 階頻率1.1163 Hz）逼近固有頻率。總耗時10.35 s，其中使用迭代法計算2794 個頻率點耗時7.20 s，使用稀疏LU 分解法計算6 個頻率點耗時3.15 s。微振動的計算成本得到了很好的控制，使未來針對微振動響應的優化提供了保障。

迭代法的迭代次數如圖2 所示。可以看出，在小于40 Hz 以下的中低頻段，迭代次數較高；在大于40 Hz 以上的中高頻段，迭代次數較低。這說明了迭代初始值在中高頻段與最終迭代值相近，其物理意義為無控制系統的頻率響應與有控制系統的頻率響應幾乎相等。即從數學角度和物理角度相互印證了控制系統在大型航天器微振動抑制中對中高頻段的效果較弱的結論。

圖2 數值計算迭代次數Fig.2 Number of iterations for numerical computation

4 結語

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