類X-51A飛行器縱向機動數值虛擬飛行仿真

王勝，王強，林博希，閻超，*

（1.北京航空航天大學 航空科學與工程學院，北京100083； 2.中國航天空氣動力技術研究院，北京100074）

吸氣式高超聲速飛行器由機體和超燃沖壓發動機一體化組成，其構型往往較為復雜［1-2］。飛行器高速飛行時，經常需要快速機動到一個較大的迎角，以實現快速的爬升。在迎角快速拉升的過程中，會產生大范圍的流動分離，使得飛行器的氣動特性出現強烈的非定常特性，氣動力遲滯效應明顯，這對控制系統的設計帶來較大挑戰。傳統基于數據庫或氣動力模型的飛行仿真，割裂了氣動、運動和控制之間的耦合關系，不能準確描述飛行器機動過程中復雜的氣動特性和運動規律，以此設計出的控制律可能難以取得令人滿意的效果。因此，需要發展一種更為先進的方法以精確模擬吸氣式高超聲速飛行器的快速機動過程。

近年來，隨著計算機技術的快速發展，數值虛擬飛行技術成為模擬飛行器機動飛行的一種新選擇［3］。數值虛擬飛行技術是一種將計算流體力學（CFD）、剛體動力學（RBD）和飛行控制系統（FCS）耦合在一起的高保真計算方法，該方法考慮了流場的非定常特性以及氣動特性和運動特性的耦合效應，可以更為精確地獲得飛行器的閉環響應特性，受到越來越多的關注［4-10］。

數值虛擬飛行的關鍵技術主要有3點：一是耦合求解CFD/RBD方程，軟件需要具備處理網格運動和變形的能力；二是姿態控制律的設計；三是控制系統與CFD系統的耦合求解方法。這3個關鍵技術中，耦合求解CFD/RBD方程、處理網格運動和變形的技術在過去幾十年中得到了長足的發展，目前已經比較成熟［11-13］。姿態控制律的設計在導航控制領域也研究較多，方法相對成熟［14-15］。而控制系統與CFD系統的耦合求解方法在國內外研究比較少，已有的關于虛擬飛行技術的研究大都采用簡化的控制律設計控制器，并將設計好的控制器以CFD 代碼的形式嵌入到CFD模塊中，達到控制目的［16］。這種做法能夠實現運動/流動/控制三者之間的耦合仿真，從而獲得比傳統方法更精確的閉環響應特性。但是也存在一些不足，比如在CFD軟件中實現控制功能編程比較復雜；每次調試更改控制器都要重新編譯CFD軟件，過程繁瑣；通常只能采用簡化的控制律函數設計控制器，只能實現較為簡單的控制功能等。

為了降低控制系統與CFD系統的耦合難度，本文基于現代軟件分布式、模塊化的發展趨勢，使用在航空航天器導航與姿態控制等領域得到廣泛應用的商業軟件Simulink實現姿態控制功能，耦合自研CFD軟件MICFD［17-18］，建立了Simulink/MICFD數值虛擬飛行仿真平臺。利用該仿真平臺，對類X-51A外形吸氣式高超聲速飛行器進行了縱向機動閉環數值仿真，通過與工程仿真結果對比，研究了運動和氣動耦合情況下非定常效應對飛行器控制響應的影響。此外，利用該平臺，還進行了一些簡單的應用，研究了縱向機動過程中舵回路時間常數對控制性能的影響。

1 數值虛擬飛行仿真平臺的實現

為了建立Simulink/MICFD數值虛擬飛行仿真平臺，首先要解決2個軟件之間的數據傳輸問題，在保證Simulink和MICFD同步運行的前提下，實現兩者之間穩定高效的數據傳輸。

Simulink可以利用MATLAB的各種命令和庫函數，實現復雜的控制仿真任務，再利用MATLAB的RTW 模塊可以將生成的Simulink模型轉變為可以直接運行的C++程序，這就為Simulink與其他應用程序的耦合提供了技術途徑。

遠程過程調用（Remote Procedure Call，RPC）［19］是一種通過網絡從遠程計算機程序上請求服務，而無需了解底層網絡技術的協議。RPC采用客戶機/服務器模式，實現進程間的同步機制，為用戶提供請求/應答的通信方式。發出請求的程序稱為客戶機，提供服務的程序稱為服務器。gRPC是由Google公司主導開發的一款語言中立、平臺中立、開源的RPC框架，支持C、C++、Python、java等多種語言版本［20］。gRPC客戶端和服務端可以在多種環境中運行和交互，由于其跨平臺、跨語言的特點，因此能夠實現不同軟件系統之間的通信，如圖1所示。

圖1 不同環境下gRPC實現過程Fig.1 Schematic diagram of gRPC implementation process in different environments

根據上述分析，可以基于gRPC網絡通信原理，在Simulink和MICFD軟件中分別建立服務器和客戶端，通過遠程調用實現2個軟件之間高效穩定的數據通信，并在此基礎上建立Simulink/MICFD數值虛擬飛行仿真平臺。如圖2所示，其具體實現步驟如下：

步驟1啟動Simulink，根據需要設計控制器，建立仿真程序，并生成可獨立運行的C++程序。

步驟2啟動MICFD，讀入網格和初始計算條件，進行定常計算，得到初始的計算流場。

步驟3進行非定常求解，計算當前時刻的力和力矩。

步驟4RBD模塊通過求解飛行力學方程，計算飛行器的姿態，將當前的飛行姿態與期望的姿態進行比較，得到姿態偏差，并將偏差信息傳遞到gRPC客戶端。

步驟5gRPC客戶端通過網絡通信向gRPC服務器發送請求，將偏差信息傳遞到gRPC服務器。

步驟6gRPC服務器接收到請求之后，調用步驟1生成的C++project，根據設計的控制器，求出每個控制面的偏轉指令，然后將舵偏指令通過網絡通信傳遞到gRPC客戶端。

步驟7CFD模塊獲得gRPC客戶端得到的舵偏指令，并根據當前時刻的姿態進行網格移動和網格重疊，準備下一物理時間步的非定常計算。

步驟8重復步驟3～步驟7，直到物理時間推進結束。

圖2 Simulink/MICFD數值虛擬飛行仿真平臺流程圖Fig.2 Flowchart of Simulink/M ICFD numerical virtual flight simulation platform

2 數值虛擬飛行計算方法

2.1 流動控制方程與求解方法

流動控制方程為三維可壓縮Navier-Stokes方程，無量綱的守恒形式為

式中：Q為守恒變量；F、G和H 為對流通量；FV、GV和HV為黏性通量；ξ、η和ζ為3個貼體坐標系方向；t為無量綱時間；Re∞為自由來流雷諾數。

流場求解采用基于結構重疊網格的有限體積法。空間離散上，無黏通量使用Roe通量差分分裂格式求解，黏性通量使用二階中心差分格式進行離散；時間推進上，采用雙時間步LU-SGS（Lower-Upper Symmetric Gauss-Seidel）方法，湍流模型采用兩方程的SST模型假設。在涉及到動態網格計算時，由于離散后的方程中含有網格體積對時間的導數項，為避免網格變化引入的額外誤差，滿足幾何守恒律十分重要［21］。由于本文使用的是剛性重疊網格技術，在計算過程中不涉及網格體積的變化，并且本文采用的空間離散方式與坐標變換格式相匹配，極大地緩解了幾何守恒律的影響，因此在本文的研究中未對幾何守恒律做單獨處理。

2.2 行力學方程求解

飛行器的運動可以分解為兩部分：質心的平動和繞質心的轉動。慣性坐標系下的質心平動動力學方程組可以表示為

式中：m為飛行器的質量；V為速度矢量；Fa為作用在飛行器上的空氣動力矢量；Fe為外力矢量，如推力等；Fg為重力矢量。體軸坐標系下，繞質心轉動的動力學方程可以表示為

式中：Hc為飛行器相對于質心的動量矩矢量；ω為角速度矢量，其分量為（ωx，ωy，ωz）；M 為力矩矢量。體軸坐標系下，繞質心轉動的運動學方程可以表示為

式中：（φ，θ，γ）為飛行器3個方向的姿態角。按式（4）方程可以用龍格-庫塔方法進行積分求解，并采用雙歐法克服方程的奇異性［22］。

2.3 姿態控制系統設計

姿態控制系統的功能是根據指令自動調整飛行器的飛行姿態，使其能夠保持穩定或者按照預定指令繞其質心轉動［23］。本文所研究的縱向姿態控制的實現如圖3所示。圖中：αc為目標迎角；α為當前時刻迎角；e為姿態角誤差；δe為控制器輸出的舵偏指令；δc為舵機輸出的舵偏角。圖3（a）為工程上常用的基于靜態數據庫的方法，該方法根據已有的靜態數據庫通過插值的方法獲得不同飛行姿態下的氣動力數據，以此氣動力數據進行飛行力學求解以獲得新的飛行姿態。圖3（b）為本文所采用的耦合CFD和控制系統的一體化模擬方法，該方法通過求解飛行器運動過程中的實時流場以獲取更為精確的非定常氣動力來代替工程方法中插值得到的氣動力，由于考慮了流場的非定常特性，其計算結果更為真實、可靠。

理論上，在數值虛擬飛行技術中可以使用任意的控制器以滿足不同的控制需求，復雜控制器的實現以及性能考察不是本文的研究重點。因此，簡單起見，在本文的研究中，控制器采用常用的PID控制器，其輸入輸出之間的關系可以表示為

圖3 縱向姿態控制Fig.3 Longitudinal attitude control

式中：e（t）為t時刻的誤差；Kp、Ki和Kd分別為比例系數、積分系數和微分系數。

舵回路采用基于位置反饋的硬反饋式回路，其近似模型可用一階慣性環節表示［24］。

式中：K為反饋增益，一般取值為1；δc（s）為δc（t）在變換后空間的表述，δc（t）為t時刻控制舵偏指令值；Td為舵回路的時間常數，反應了舵機的遲滯特性，其值的大小受舵機的功率等自身因素影響。

3 算例驗證

為了考核MICFD軟件中CFD/RBD耦合計算能力，選取三維機翼掛載分離模型，對程序進行驗證。該模型由3部分組成，即機翼、掛架以及外掛物，其幾何外形具體參數詳見文獻［25］。計算條件為：Ma=0.95、H=7.92 km、α=0°。

圖4和圖5分別為分離過程中掛載物的氣動力、力矩以及質心位移和角位移隨時間的變化曲線。可以看出，各曲線的計算結果與實驗結果［25］吻合良好，驗證了本文CFD/RBD耦合計算方法的準確性，表明軟件具備良好的非定常多體相對運動數值模擬能力。

圖4 掛載物氣動力系數、氣動力矩系數隨時間變化曲線Fig.4 Time history of aerodynam ic coefficients and aerodynamic moment coefficients for store

圖5 掛載物線位移、角位移隨時間變化曲線Fig.5 Time history of linear and angular displacement for store

4 類X-51A飛行器縱向機動過程閉環模擬

4.1 計算模型和網格

本文首先基于X-51A飛行器的相關資料［1，26］，獲得其幾何外形的主要參數，采用反向建模技術建立了一種與X-51A外形相似的高超聲速飛行器模型，如圖6所示。本文主要研究飛行器縱向機動過程，模型左右對稱且不考慮側滑，為節省計算量采用半模計算。網格生成時，飛行器本體網格和控制舵面網格分別獨立生成，并采用動態網格重疊的方式實現兩者的相對運動。網格總量約1 600萬。機動過程中不同時刻控制舵與機身的動態重疊邊界示意圖如圖7所示。

圖6 計算模型及對稱面網格Fig.6 Computationalmodel and symmetry plane grids

圖7 不同時刻控制舵與機身動態重疊邊界示意圖Fig.7 Schematic diagram of dynamic overlapping boundary of control rudder and airframe at differentmoments

4.2 靜態氣動特性

選取的計算狀態為典型的進氣道起動狀態。來流馬赫數為4.8，來流單位雷諾數為1.2×107。力矩計算時，為方便起見，參考長度取Lref=1m，參考面積取Sref=1m2。坐標原點設置在飛行器前緣縱向對稱面處。力矩參考點（x，y，z）=（1.8，-0.5，0）m。

圖8顯示了不同迎角和舵偏角下飛行器的俯仰力矩，其中點畫線代表的是全機俯仰力矩，實線是控制舵產生的俯仰力矩。可以看出，在所計算的迎角范圍內，除20°舵偏角之外，其余舵偏下俯仰力矩的線性度比較好。而全機俯仰力矩具有較明顯的非線性。圖9給出了各迎角和舵偏角下除控制舵以外的機身部分的俯仰力矩。可以看出，與傳統的軸對稱外形或者升力式外形飛行器不同，該外形機身部分的俯仰力矩呈現明顯的非線性，說明由于進氣道和內流道的存在，迎角的改變對氣動特性的影響更加復雜。

圖8 全機及控制舵俯仰力矩Fig.8 Pitching moment of whole aircraft and control rudder

圖9 除控制舵以外機身的俯仰力矩Fig.9 Pitching moment of airframe except control rudder

4.3 縱向機動過程模擬

本節對建立的類X-51A模型進行縱向機動仿真。計算的來流條件與靜態計算時相同，初始的來流迎角為0°，控制指令為：以階躍響應的方式機動到10°迎角，并保持穩定。分別采用圖3（a）介紹的工程方法和圖3（b）介紹的耦合方法進行計算。由4.2節計算可知，初始狀態下控制舵面的配平偏轉角為-11°，閉環模擬時以該舵偏角作為初始舵偏，即假定開始機動時，飛行器處于平衡狀態（后文中的舵偏角均指在此基礎上的相對舵偏角）。

通過試湊，選取一組參數：Kp=2，Ki=1.5，Kd=0.6作為控制器的增益。舵回路中，反饋增益K=1，舵回路時間常數Td=0.1。采用工程方法計算時，飛行器不同迎角和舵偏角下的氣動力由靜態數據庫通過插值得到。非定常計算時，物理時間步長設為1ms，子迭代殘差指標為0.01，同時限定最大子迭代步數為50步，子迭代CFL數為3.0。

圖10給出了采用耦合方法模擬時不同時刻飛行器周圍的流場結構，包括壁面壓力分布、軸向不同截面的馬赫數等值線以及控制舵附近的流線等。可以看出，相比于圖10（a）的初始定常狀態，飛行器在拉起過程中迎風面和背風面出現明顯差異，控制舵附近流線偏轉，出現了較強的非定常特征。

圖10 機動過程中典型時刻的流場結構（T d=0.1）Fig.10 Typical-moment flow field structure during maneuvering process（T d=0.1）

圖11～圖13分別表示了2種方法得到的迎角、舵偏角以及俯仰力矩的響應過程。可知，2種方法均可以實現控制指令要求的機動過程，最終的配平舵偏角均保持在-5.6°附近，但是兩者之間也存在一些差異。從迎角響應過程來看，耦合方法計算得到的超調量（16%）略高于工程方法（13.4%），耦合方法計算得到的調節時間為2.01 s，遠高于工程方法得到的1.48 s，說明工程算法給出的結果可能低估了控制系統的遲滯特性。從舵偏角響應過程來看，盡管最終的配平舵偏相同，但是兩者在峰值位置存在較大差異，這正反映了傳統的基于靜態氣動力數據庫的飛行仿真在模擬快速機動過程中的缺陷，體現了流場非定常效應對控制過程的影響，這一點通過圖13中俯仰力矩在峰值處的差異也可以看出。因此，對于類似X-51A、外形復雜的高超聲速飛行器，由于非定常效應明顯，很有必要采用基于耦合方法的數值虛擬飛行技術來研究飛行器的機動過程、評估飛行控制律。

圖11 迎角響應過程Fig.11 Response process of angle of attack

圖12 舵偏角響應過程Fig.12 Response process of rudder deflection angle

圖13 俯仰力矩時間歷程Fig.13 Time history of pitching moment

4.4 舵回路時間常數的影響

舵回路系統也稱伺服系統，是飛行控制系統中的重要組成部分。舵機響應的遲滯特性對飛行器的控制系統有著顯著影響。因此，對控制系統在舵機不同響應特性下的控制性能進行評估和分析，為控制律的設計提供參考具有重要的意義。在本節中，通過改變舵回路系統的時間常數Td進行數值虛擬飛行仿真，以評估其對飛行操縱過程的影響。其他仿真條件與4.3節一致。

圖14～圖16給出了迎角、舵偏角以及俯仰力矩在不同舵回路時間常數下的響應過程。可以看出，不考慮舵回路（即Td=0）時，舵偏操縱是理想的，即δc=δe。控制舵一開始就存在一個較大的偏轉角（約-20°）；而在考慮舵回路時，由于舵機的遲滯特性，舵偏角是從0°開始逐漸增大的，這也導致2種情況下飛行器的機動過程出現較大差異。表1列出了不同時間常數下控制系統的性能指標，包括系統的延遲時間（td）、上升時間（tr）、峰值時間（tp）、調節時間（ts）以及超調量（σ%）。可以看出，隨著Td的增加，系統的調節時間和超調量顯著增加，動態性能變差。值得注意的是，當Td=0.2時，超調量達到31.6%，并且縱向姿態在數值模擬的時間歷程中未能收斂到指定的穩定姿態。可以預見，繼續增大Td可能會使系統難以收斂甚至發散。

從以上分析可以看出，舵回路的時間常數對控制過程有重要影響。在滿足舵機功率限制的前提下，通過設計適當的舵回路來減小Td，可以降低超調量，提高響應速度，改善系統的動態性能。

圖14 不同時間常數下迎角響應過程Fig.14 Response process of angle of attack under different time constants

圖15 不同時間常數下舵偏角響應過程Fig.15 Response process of rudder deflection angle under different time constants

圖16 不同時間常數下俯仰力矩時間歷程Fig.16 Time history of pitching moment under different time constants

表1 不同T d 下控制系統的性能指標Tab le 1 Perform ance indexes of control system under d ifferent T d

5 結 論

本文從現代軟件分布式、模塊化的發展趨勢出發，基于 gRPC 網絡通信建立了 Simulink/MICFD數值虛擬飛行仿真平臺，有效降低了控制模塊和CFD模塊的耦合難度。利用該仿真平臺，對類X-51A外形的吸氣式高超聲速飛行器進行了縱向機動閉環數值仿真。得到如下結論：

1）本文采用的遠程過程調用的方法在形式上完全分割了控制和CFD這2個模塊，但是在邏輯上又將2者緊密耦合在一起。相比代碼級耦合，基于此方法建立的數值虛擬飛行仿真平臺能夠充分發揮控制和CFD軟件各自的優點，顯著降低了多學科耦合的難度。

2）仿真結果表明，該數值平臺具備針對復雜飛行力學行為和控制響應特性的先進數值模擬能力，能夠實現對飛行器有控機動過程的精細化模擬。

3）對于類X-51A外形的吸氣式高超聲速飛行器，在縱向拉起時，工程算法給出的結果可能不能完全反映非定常效應的影響。有必要采用基于耦合方法的數值虛擬飛行技術來研究飛行器的機動過程、評估飛行控制律。

4）舵回路的時間常數對控制系統的性能有重要影響，減小時間常數可以降低超調量，提高響應速度，改善系統的動態性能。