飛行器運動控制建模及紅外輻射圖像仿真

劉勇求，金忠慶

飛行器運動控制建模及紅外輻射圖像仿真

劉勇求1，金忠慶2

（1.廣東科技學院，東莞 523083；2. 空軍航空大學，吉林 長春 130022）

實時生成飛行器在飛行過程中的紅外圖像對于研究飛行器的紅外特性具有重要意義。飛行器在飛行過程中，由于其飛行速度、高度、姿態(tài)等參數(shù)變化性較大，且機體本身容易受到環(huán)境、背景等輻射的影響，因此較難準確地仿真機體的溫度場分布。由于飛行器機體溫度場分布與飛機的運動狀態(tài)密切相關(guān)，因此文中首先建立了決定其姿態(tài)和速度的控制模型，并對飛行器面源熱平衡方程進行建模研究。然后，應(yīng)用CFD（Computational Fluid Dynamics）對飛行器的絕熱壁溫度場進行求解，并給出了機體絕熱壁溫度的高度簡化算法。最后仿真了不同飛行器在不同波段、不同高度、不同速度下的紅外圖像。

紅外圖像；面源熱平衡方程；CFD；絕熱壁溫度

0 引言

飛行器在高速飛行過程中其機體蒙皮表面與周圍的大氣相互摩擦，從而產(chǎn)生大量熱量，同時發(fā)動機產(chǎn)生的熱量也會通過熱傳導(dǎo)和熱輻射傳遞到機體蒙皮。因此在飛行過程中，飛行器機體的溫度明顯高于周圍大氣，其紅外特征較為明顯。機體蒙皮各面源的溫度差異較大，其機翼前緣、機頭、座艙蓋上部由于與空氣摩擦劇烈溫度較高，機身尾部貼近發(fā)動機噴口處由于受發(fā)動機熱效應(yīng)的影響其蒙皮溫度也較高，機體其他處溫度較低，因此各面源紅外圖像的灰度分布也不一樣。由于飛行器在飛行過程中其速度、高度、姿態(tài)變化較大，且其紅外輻射易于受到環(huán)境、背景等輻射的影響，因此較難準確地仿真機體的溫度場分布。目前，國內(nèi)外研究飛行器點源紅外特征的較多，研究機體紅外輻射特性的較少。劉娟[1]研究了飛機的紅外輻射特性；李建勛[2]將機體的紅外輻射分為機身前部和尾部輻射兩部分，分別研究其紅外特性；黃鶴松[3]研究了自燃活性金屬的熱平衡方程，并仿真了其燃燒過程中的紅外圖像。

由于飛行器在飛行過程中其飛行姿態(tài)、速度、高度以及各舵面的操縱量不斷變化，因此本文首先建立飛行器的運動與控制模型，得到其飛行過程中的各種運動狀態(tài)；然后建立飛行器面源熱平衡方程，并應(yīng)用計算流體力學（Computational Fluid Dynamics，CFD）對飛行器的絕熱壁溫度場進行了求解；最后對飛行器絕熱壁溫度場的計算結(jié)果進行高度簡化，并仿真了飛行器在不同高度、速度、波段下的紅外輻射圖像。

由于飛行器發(fā)動機紅外輻射特性方面的研究成果較多，因此本文采用文獻[4-9]中的模型，不再對其進行單獨建模。

1 飛行器運動與控制模型

飛行器紅外輻射圖像的仿真需要得到飛行器速度、姿態(tài)、舵偏角等數(shù)據(jù)，因此對飛行器的運動與控制規(guī)律進行建模十分必要。飛行器在飛行中主要受到的外力有：氣動力、自身重力以及發(fā)動機推力，這些外力通常按不同的坐標系給出，同時飛行器的運動也需要在相應(yīng)的坐標系內(nèi)進行描述，各坐標系之間的關(guān)系詳見參考文獻[10]。

1）飛行器質(zhì)心運動動力學方程

式中：V、V、V、、、分別為飛行器三軸速度、角速度分量；為發(fā)動機的推力；g為自身重力；、、、、分別為迎角、側(cè)滑角、俯仰角、滾轉(zhuǎn)角以及發(fā)動機安裝角；、、分別為氣動阻力、氣動升力以及側(cè)力；為時間。

2）飛行器繞質(zhì)心轉(zhuǎn)動動力學方程

式中：?M、?M、?M為飛行器的合外力矩；I、I、I為飛行器的三軸慣性矩；I、I、I為飛行器的慣性積，其中I＝0，I＝0。

3）飛行器繞質(zhì)心轉(zhuǎn)動運動學方程

聯(lián)立公式(1)～(3)，在知道了飛行器的空氣動力特性、發(fā)動機推力特性和發(fā)動機安裝角、初始飛行狀態(tài)和各種姿態(tài)角后，通過數(shù)值積分，即可解得飛行器在運動過程中的速度和姿態(tài)。

2 飛行器機體熱平衡方程

飛行器在高速飛行過程中機體與空氣摩擦，機體蒙皮溫度上升，其輻射亮度值主要取決于機體蒙皮表面的溫度場及機體蒙皮發(fā)射率。機體蒙皮發(fā)射率一般可當作固定值處理，因此求解機體的紅外輻射亮度關(guān)鍵在于準確地求解機體表面的溫度場。

機體的表面溫度與機體表面換熱量有以下關(guān)系：

式中：D為單位時間內(nèi)單位面積的機體蒙皮上的熱量變化值；為機體蒙皮溫度；為比熱容，則有：

D＝in＋rad＋flow＋env(5)

式中：in為蒙皮之間熱傳導(dǎo)、熱輻射產(chǎn)生的能量；rad為蒙皮向外界輻射的能量；flow為蒙皮與周圍大氣之間的換熱量；env為蒙皮接受到周圍環(huán)境輻射的能量。

2.1 蒙皮之間熱傳導(dǎo)、熱輻射產(chǎn)生的能量

蒙皮之間熱傳導(dǎo)、熱輻射的能量主要發(fā)生在發(fā)動機艙，飛行器在飛行過程中，由于發(fā)動機燃燒室內(nèi)猛烈持續(xù)的燃燒反應(yīng)，發(fā)動機的各個部件都會產(chǎn)生大量的熱量，這些熱量大部分通過發(fā)動機的冷卻系統(tǒng)進行冷卻，少部分熱量通過熱傳導(dǎo)和熱輻射的方式傳遞到相鄰的其他蒙皮上。蒙皮之間熱傳導(dǎo)、熱輻射的能量可以由下式得到：

式中：in為發(fā)動機的簡化溫度；in為發(fā)動機的發(fā)射率；sur為蒙皮表面的吸收率；in-sur為輻射角系數(shù)；c為熱傳導(dǎo)系數(shù)；in-sur為傳導(dǎo)距離；sur為蒙皮表面溫度；為斯忒藩-玻爾茲曼常數(shù)，其值為5.67×10－8W·m－2·K－4。

2.2 機體表面向外輻射熱量

由玻爾茲曼定律，黑體的輻射出射度為：

＝4(7)

將機體蒙皮簡化為發(fā)射率sur的灰體，則單位時間內(nèi)機體蒙皮自身輻射的熱量為：

rad＝－sursur4(8)

2.3 環(huán)境輻射熱量

環(huán)境輻射主要包括大氣輻射、海洋輻射、地面輻射、云層輻射、太陽輻射等，由于飛行器在高空中飛行，地面輻射對其影響不大，本文只考慮對機體蒙皮影響最大的大氣輻射和太陽輻射[11-12]。環(huán)境輻射可由下式得到：

2.4 機體蒙皮與周圍大氣的換熱量

飛行器在空中飛行過程中，其機體蒙皮表面與周圍大氣摩擦，摩擦作用使周圍氣體溫度上升。同時，在飛行過程中，飛行器蒙皮表面壓縮前方氣體，使得氣體溫度進一步升高。在飛行器的飛行過程中，機體的附面層氣體溫度較高，則機體蒙皮與附面層氣體的換熱量為：

flow＝flow(flow－sur) (10)

式中：flow為附面層氣體溫度；flow為換熱系數(shù)。換熱系數(shù)flow主要由空氣導(dǎo)熱系數(shù)a；機體特征長度以及雷諾數(shù)Re確定，其值為：

式中：Re為雷諾數(shù)；Pr為普朗特數(shù)，空氣的Pr數(shù)一般取0.72。附面層氣體的溫度flow既受飛行器外形的影響，又與飛行器的飛行速度、高度、大氣溫度、大氣密度等有關(guān)，無法精確獲得。文中采用CFD計算飛行器機體蒙皮表面的絕熱壁溫度，然后用絕熱壁溫度來近似代表flow，求解過程中設(shè)置飛行器的機體表面為絕熱壁[13]。

3 機體絕熱壁溫度計算

3.1 網(wǎng)格刨分

首先建立飛行器的三維幾何模型。為了保證模型的精確性，本文采用CATIA三維建模軟件對飛行器的三維模型進行建模，以F-22為例，建立的三維幾何模型如圖1所示。

圖1 F-22三維模型圖

利用ICEM軟件對計算域及機體表面進行網(wǎng)格劃分。本文以F-22為例，設(shè)計算區(qū)域為210m×70m×70m的長方體，飛行器中心與計算域的中心相重合，如圖2所示。邊界層網(wǎng)格第一層厚度為10－3m，邊界層一共12層，采用四面體網(wǎng)格，如圖3～4所示。

圖2 計算域示意圖

圖3 F-22機體表面網(wǎng)格圖

圖4 計算域網(wǎng)格圖

3.2 CFD流場計算

控制方程選用Navier-Stokes方程。采用有限體積的空間離散算法和空間二階精度的線性插值算法，時間離散采用二階精度的后向差分算法，計算模型選用SST湍流模型[10]。計算域外邊界設(shè)置為壓力遠場邊界，機體表面設(shè)置為無滑移絕熱壁面，進氣道入口設(shè)置為質(zhì)量入口，質(zhì)量流量由發(fā)動機氣體流量數(shù)據(jù)庫決定，尾噴口設(shè)置為壓力出口，認為噴口氣體處于完全膨脹狀態(tài)，則噴口靜壓與外界大氣壓力相同。

設(shè)置的初始條件為高度11km，壓力＝22700Pa，大氣溫度＝216.774K，飛行器飛行馬赫數(shù)Ma＝1.6。則F-22飛機絕熱壁溫度分布計算結(jié)果如圖5所示。

圖5 F-22飛機絕熱壁溫度分布圖

飛行器在穩(wěn)定飛行過程中，機體蒙皮各面源溫度不變，則公式(5)中的D＝0。因此，在得到飛行器各面源的絕熱壁溫度flow后，聯(lián)立公式(5)～(11)，便可以求解出飛行器機體蒙皮的表面溫度sur。

3.3 高度簡化算法

由于飛行器的飛行高度通常在0～20km內(nèi)，飛行馬赫數(shù)在0～2范圍內(nèi)變化，飛行過程中由于舵面的操縱，迎角和側(cè)滑角也都在變化。因此要計算飛行器在運動過程中的紅外輻射圖像，需要對飛行器在各種條件下的絕熱壁溫度進行計算，計算量十分龐大，需要對CFD計算結(jié)果進行簡化。

當飛行器飛行速度增加時，其機翼表面上的雷諾數(shù)也跟著增加，而當雷諾數(shù)增加到一定數(shù)值后，機體表面的附面層就會出現(xiàn)轉(zhuǎn)戾現(xiàn)象，從而出現(xiàn)紊流，使得機體表面的溫度分布差異很大。因此飛行器蒙皮的絕熱壁溫度分布隨其飛行速度的變化關(guān)系較為復(fù)雜、無章，無法對其進行簡化。而飛行高度主要決定著周圍大氣的環(huán)境溫度，對飛行器機體的溫度分布影響一般較小，因此本節(jié)研究飛行器機體絕熱壁溫度的高度簡化算法。

飛行器機體的駐點溫度可由下式得到：

式中：0為大氣環(huán)境溫度；為恢復(fù)溫度系數(shù)（層流取值一般為0.82，湍流為0.87）；為大氣絕熱指數(shù)，＝1.4；b為駐點溫度。設(shè)機體蒙皮的絕熱溫度與駐點溫度比為，即：

flow＝b(13)

式中：表示機體蒙皮各面源與駐點溫度的比值，也表示了機體的溫度分布情況。

設(shè)F-22飛機飛行馬赫數(shù)為0.6，則圖6為其在不同高度下機體絕熱壁溫度場的CFD計算結(jié)果。

由圖6可知，當飛行器的飛行馬赫數(shù)相同時，不同高度下飛行器的絕熱壁溫度分布基本不變。設(shè)蒙皮面源在km高度時絕熱壁溫度為T,i，T為飛行器在km高度時的駐點溫度，則蒙皮面源的簡化結(jié)果為：

式中：T,i￠為面源的絕熱壁溫度簡化計算結(jié)果；,i為采用簡化算法后的面源的絕對誤差。采用簡化算法得到的誤差對比分布圖如圖7和圖8所示。

圖7 各高度下eh,i的分布圖（Ma＝0.6）

圖8 各高度下eh,i的分布圖（Ma＝1.5）

由圖7和圖8可知，當飛行馬赫數(shù)為0.6時，絕對誤差,i的值全部小于1K；當飛行馬赫數(shù)為1.5時，絕對誤差,i的值全部小于4.5K，且只有當高度為3km時，少部分面源的絕對誤差,i值大于3K，其余條件下全部小于1.8K。由圖9可知，飛行高度在10km附近誤差最小，飛行高度在海平面和20km時，誤差最大，但整體誤差仍較小。因此采用公式(14)的簡化算法計算得到的絕熱壁溫度結(jié)果較為準確，能夠滿足計算精度的要求。

4 機體紅外輻射圖像仿真

令D＝0，則由公式(5)及飛行器的絕熱壁溫度CFD計算結(jié)果，可以求得飛行器蒙皮任意面源的溫度，從而得到飛行器在飛行過程中的溫度分布結(jié)果。本節(jié)應(yīng)用VC＋＋語言，借助OpenGL繪圖工具繪制飛行器的紅外輻射圖像。

設(shè)F-22飛機水平飛行，飛行馬赫數(shù)為0.6。飛行器平尾舵偏角＝5°，垂尾和副翼的舵偏角＝0°、＝0°，則其在不同高度下的紅外圖像如圖10所示。

設(shè)F-22飛機水平飛行，飛行馬赫數(shù)為1.0。飛行器平尾舵偏角＝3°，垂尾和副翼的舵偏角＝0°、＝0°，則其在不同高度下的紅外圖像如圖11所示。

設(shè)F-16飛機水平飛行，飛行馬赫數(shù)為0.6時，平尾舵偏角＝7°，垂尾和副翼的舵偏角＝0°、＝0°；飛行馬赫數(shù)為1.0時，平尾舵偏角＝4°，垂尾和副翼的舵偏角＝0°、＝0°，則其在不同高度下的紅外圖像仿真結(jié)果如圖12和圖13所示。

由圖10～圖13的仿真結(jié)果可知，飛行器紅外圖像的亮度隨著飛行高度的上升而下降，隨著飛行速度的增加而增加，同時8～12mm波段的圖像亮度明顯高于3～5mm波段。亮度越高的點表示此處的紅外輻射亮度越強，相應(yīng)的溫度也越高。對比相同條件下的仿真結(jié)果可知，F(xiàn)-22飛機的紅外輻射強度略大于F－16，這是由于文中仿真未考慮F-22飛機隱身涂層，且F-22飛機機身蒙皮面積大于F-16所致。

圖9 絕對誤差最大值隨高度變化的曲線

圖11 F-22紅外輻射圖像（Ma＝1.0）

圖12 F-16紅外輻射圖像（Ma＝0.6）

5 結(jié)論

本文對飛行器的紅外圖像進行仿真建模，由于飛行器在運動過程中其姿態(tài)的變化主要受控制模型的影響，因此文中首先建立了飛行器運動與控制模型；然后研究了飛行器面源熱平衡方程，并應(yīng)用CFD對飛行器的絕熱壁溫度場進行了求解；最后研究了飛行器絕熱壁溫度的高度簡化算法，并對飛行器的紅外輻射圖像進行了仿真計算。所得到的主要結(jié)論有：

1）飛行器蒙皮的絕熱壁溫度場可以進行高度簡化處理，簡化后大部分面源的絕熱壁溫度與CFD計算結(jié)果誤差在3K以內(nèi)。

2）飛行器飛行高度在20km以下時，其紅外圖像的亮度隨著飛行高度的上升而下降，隨著飛行速度的增加而增加。

3）8～12mm波段的圖像亮度明顯高于3～5mm波段。

[1] 劉娟, 龔光紅, 韓亮, 等.飛機紅外輻射特性建模與仿真[J]. 紅外與激光工程, 2011, 40(7): 1219-1213.

LIU Juan, GONG Guanghong, HAN Liang, et al. Modeling and simulation of airplane infrared characteristic[J]., 2011, 40(7): 1219-1213.

[2] 李建勛, 童中翔, 王超哲, 等. 飛機紅外輻射特性分析與易損性評估[J]. 激光與紅外, 2013, 43(2): 180-185.

LI Jianxun, TONG Zhongxiang, WANG Chaozhe, et al. IR radiation characteristics analysis and vulnerability assessment of aircraft[J]., 2013, 43(2): 180-185.

[3] HUANG Hesong, TONG Zhongxiang, WANG Chaozhe, et al. Numerical and Experimental study of pyrophoric activated metal Mg surface combustion characteristics[J]., 2018, 5: 172064.

[4] 汪中賢, 史建軍, 張正武. 噴氣式飛機尾焰紅外輻射的理論計算與仿真[J]. 火力與指揮控制, 2011, 36(11): 83-86.

WANG Zhaongxian, SHI Jianjun, ZHANG Zhangwu. Theoretical Calculation and Simulation of Infrared Radiation of Jet Aircraft Plume[J]., 2011, 36(11): 83-86.

[5] 陳衛(wèi), 汪中賢, 馬東輝, 等. 非均勻熱氣體紅外輻射特性計算與仿真[J]. 紅外與激光工程, 2010, 39(1): 17-21.

CHEN Wei, WANG Zhongxian, MA Donghui, et al. Calculation and simulation of infrared radiation characteristics of non-uniform hot gas[J]., 2010, 39(1): 17-21.

[6] 高思莉, 湯心溢. 空中飛行目標尾焰紅外輻射信號的建模與仿真[J]. 光電工程, 2007, 34(8): 25-27.

GAO Sili, TANG Xinyi. Building model of the plume released from the flying machine and simulation [J]., 2007, 34(8): 25-27.

[7] 范仁鈺. 下遮擋二元噴口的氣動和紅外輻射特性的數(shù)值研究[D]. 南京: 南京航空航天大學, 2009.

FAN Renyu. Numerical Investigation on the Aerodynamics and Infrared Radiation Characteristics of 2-Dimensional Nozzles with Sheltering Baffles[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009.

[8] 李建勛, 童中翔, 王超哲, 等. 發(fā)動機熱噴流紅外輻射計算與仿真[J]. 光譜學與光譜分析, 2013, 33(1): 7-13.

LI Jianxun, TONG Zhongxinag, WANG Chaozhe, et al. Calculation and Simulation on Infrared Radiation of Hot Jet from Engine[J]., 2013, 33(1): 7-13．

[9] HU Chongwei, Joe Der, Jr, et al. A simple 2D-nozzle plume model for IR analysis[R]. AIAA-80-1808, 1980.

[10] HUANG Hesong, TONG Zhongxiang, CHAI Shijie, et al. Experimental and numerical study of chaff cloud kinetic performance under impact of high-speed airflow[J]., 2018, 31(11): 2080-2092.

[11] 張志波, 童中翔, 王超哲, 等. Lowtran在紅外成像仿真中的應(yīng)用[J]. 計算機應(yīng)用與軟件, 2013, 30(7): 24-26.

ZHANG Zhibo, TONG Zhongxiang, WANG Chaozhe, et al. The Application of Lowtran in Infrared Imaging Simulation[J]., 2013, 30(7): 24-26．

[12] 吳北嬰, 李衛(wèi), 陳洪濱, 等. 大氣輻射傳輸實用算法[M]. 北京: 氣象出版社, 1998.

WU Beiying, LI Wei, CHEN Hongbin, et al.[M]. Beijing: Meteorological Publishing House, 1998.

[13] 李劉苗. 海天背景目標紅外輻射特性研究[D]. 西安: 西安電子科技大學, 2011.

LI Liumiao. Research on Infrared Radiation Characteristic of Object Under Sky and Sea Background[D]. Xi’an: Xidian University, 2011.

Simulation of Air Vehicle Motion Control Models and Infrared Images

LIU Yongqiu1，JIN Zhongqing2

(1.,523083,; 2.,130000,)

The generation of infrared (IR) images of air vehicles in real time during flight is very important for research on the IR characteristics of air vehicles. The velocity, altitude, and attitude of air vehicles vary during flight. Moreover, the air vehicle body is affected by environmental and background radiation. Hence, it is very difficult to obtain an accurate temperature distribution of the air vehicle body. Because the temperature distribution is related to the state of motion, control models that determine the attitude and velocity are established, along with the surface heat balance equations. The adiabatic wall temperature distribution is solved by computational fluid dynamics, and simplified altitude algorithms are obtained. Finally, IR images of the air vehicle in different bands and at different altitudes and velocities are simulated.

IR image, surface heat balance equations, CFD, adiabatic wall temperature

TJ7；TP391

A

1001-8891(2020)09-0863-10

2019-12-31；

2020-06-01.

劉勇求（1980-），女，湖南湘潭人，碩士，副教授，研究方向為智能控制。

廣東省普通高校特色創(chuàng)新類項目（2018KTSCX262）。