基于新型智能探頭的大氣數據系統算法研究

劉 駿，鄭 偉，王 禹

(航空工業沈陽飛機設計研究所，遼寧 沈陽 110035)

大氣數據系統是一種具有綜合性和高精度特點的大氣數據測量解算系統[1-3]，可以提供高度Hp、空速Vc、馬赫數Ma、迎角αt和側滑角βt等飛行器飛行所需的重要飛行參數，對于保障飛行器的飛行安全起著至關重要的作用。隨著我國科學技術的不斷發展進步，相關的技術工藝水平也得到了進一步的提高[4]，大氣數據系統已由最早期的多個分立式傳感器組合，發展到由多個智能探頭組成的分布式大氣數據系統[5]。

在基于智能探頭的分布式大氣數據系統中，每一個或每一組智能探頭都相當于一套獨立的大氣數據測量輸出系統，可以根據飛機測量余度的需要，在機身不同位置處分別安裝數只，組成一套分布式大氣數據系統。根據測量形式的不同，目前常用的智能探頭可分為3種[6]，分別是風標式智能探頭、空速管式智能探頭和錐形智能探頭，如圖1所示。

本文主要設計了一種基于半圓柱形智能探頭(下文簡稱為半圓柱探頭)的分布式大氣數據系統的算法模型(該種半圓柱探頭由測壓傳感器、智能電路和數據接口等部分組成)，并開展了基于此算法模型的仿真驗證工作，為后續大氣數據系統的設計研究工作提供一種新的思路。

1 半圓柱探頭外形及安裝信息

該種半圓柱探頭整體呈半圓柱形，圓柱端為安裝端，自由端為半圓柱型，外側靠近自由端處開4個測壓孔，其中一個測壓孔在半圓柱體臺階平面上，另外3個測壓孔均勻分布在背部曲面上，半圓柱探頭外形示意圖如圖2所示。

圖2 半圓柱探頭布局示意圖

半圓柱探頭在機上安裝后，通過CFD(Computational Fluid Dynamics，計算流體力學)軟件計算得到機身附近局部流場信息，選擇半圓柱探頭安裝位置[7-8]，將每兩個半圓柱探頭為一組，垂直蒙皮對稱安裝在飛行器前機身左下和右下腹部處，探頭圓柱面正對來流方向，臺階平面處于背風面。

2 半圓柱探頭的大氣數據系統算法建模

2.1 算法模型基本原理

分布式大氣數據系統測量大氣數據的基本原理[9-10]是：通過風洞試驗和CFD計算，標定得到在不同狀態下智能探頭感受到的局部來流壓力、角度和飛機遠端自由來流的壓力、角度之間的定量關系，之后再由大氣解算部件利用局部壓力、角度計算得到表征自由流的大氣參數。按照基于半圓柱探頭的大氣數據系統方案構型[11]，系統應包括對稱安裝于前機身左右表面的4個半圓柱探頭(兩兩對稱)，每個探頭(橫截面外形如圖3所示)感受并測量輸出4個壓力給大氣解算部件。

圖3 半圓柱探頭橫截面示意圖

具體來說，機頭前遠方自由來流Pt(QC)、PS、Ma、αt、βt受機身擾動，氣流流經左、右探頭所在位置處時，其相對左、右探頭的局部來流分別為PLt(QLC)、PLS、MaL、θL和PRt(QRC)、PRS、MaR、θR。左、右局部來流再流經左、右探頭時，左、右探頭上各4個測壓孔感受的壓力分別為PL1、PL2、PL3、PL4和PR1、PR2、PR3、PR4。

大氣數據系統算法軟件模塊駐留在大氣解算部件的每個通道中，以其中一個通道為例，其接收對稱位置兩個探頭的8個壓力信號為PL1、PL2、PL3、PL4、PR1、PR2、PR3、PR4，先基于半圓柱體探頭本體的氣動模型，利用左、右半圓柱探頭各自4個點壓力解算得到各自所在位置機身局部流場的QLC、PLS、MaL、θL和QRC、PRS、MaR、θR。

再根據探頭所在位置局部流場的氣動模型，利用左、右探頭局部來流角均值求解真迎角αt，利用左、右探頭局部來流角差值求解真側滑角βt，取左、右探頭局部流場動壓均值為指示動壓QCi，取左、右探頭局部流場靜壓均值為指示靜壓PSi，利用QCi和PSi之比求解指示Mai。再以αt、βt和Mai為基礎，根據CFD計算獲取的半圓柱探頭所在位置局部流場的氣動模型，完成由指示值到真實值的迭代計算，得到當前飛行狀態下的QC、PS、Ma、αt和βt。

2.2 算法模型相關參數定義

本節給出氣動建模所需部分樣本數據涉及到的各參數定義。

① 左、右探頭上各測壓點的壓力系數定義：

CpLi(MaL,θL)=(PLi-PLS)/QLC，i=1，2，3，4

(1)

CpRi(MaR,θR)=(PRi-PRS)/QRC，i=1，2，3，4

(2)

② 左、右探頭位置處機身局部流場靜壓誤差壓力系數：

CpLS(Ma,αt,βt)=(PLS-PS)/QC

(3)

CpRS(Ma,αt,βt)=(PRS-PS)/QC

(4)

③ 左、右探頭位置處機身局部流場總壓誤差壓力系數：

CpLt(Ma,αt,βt)=(PLt-Pt)/QC

(5)

CpRt(Ma,αt,βt)=(PRt-Pt)/QC

(6)

④ 左、右機身局部流場來流角與馬赫數Ma、真迎角αt、真側滑角βt關系：

θL=θL(Ma,αt,βt)

(7)

θR=θR(Ma,αt,βt)

(8)

2.3 氣動模型設計與分析

從2.1節的內容中可知，大氣數據系統算法建模的核心分為兩方面，一是獲取半圓柱探頭本體4個測壓點壓力系數樣本數據，并建立相應測壓點間壓力系數差與局部來流動、靜壓和角度之間關系；二是獲取半圓柱探頭所在位置處局部流場的氣動模型樣本數據，并建立局部來流和自由來流之間的映射關系[12-13]。在確定探頭外形和探頭在機上的安裝位置和安裝角度之后，上述壓力系數樣本數據即可通過風洞試驗和CFD計算獲取。

2.3.1 半圓柱探頭本體氣動模型

以飛機左側探頭為例，首先根據CFD計算數據，做出探頭上P1、P2、P3和P4四個測壓點處不同馬赫數下壓力系數隨局部來流角變化的曲線，圖像如圖4和圖5所示。

圖4 0.4馬赫數下各測壓點處壓力系數隨局部迎角變化曲線

圖5 0.8馬赫數下各測壓點壓力系數隨局部迎角變化曲線

通過對圖4和圖5中的曲線進行分析，可以發現在不同馬赫數下，P2點處壓力系數對局部流場來流角變化不太敏感，在不同馬赫數下數值雖隨之變化但基本能保持穩定，結合圖3所示流場來流和探頭處位置關系示意圖，可以確定用P2測壓點處測量得到的壓力值可以作為探頭所在處的局部靜壓的指示靜壓使用。同時，P1、P3和P4三個測壓點處的壓力系數在不同馬赫數下的隨局部來流角變化時，總有一個點的壓力系數接近1，說明此時該點近似于正對來流方向，該測壓點感受到的壓力值可以作為探頭所在處的局部總壓的指示總壓使用，所以可得到局部流場指示靜壓：

PLSi=PL2

(9)

局部流場指示動壓：

QLCi=PLmax-PLSi

(10)

式中：PLmax為PL1、PL3、PL4中的最大值。

局部流場指示馬赫數：

MaLi=f(QLCi/PLSi)

(11)

式中：MaLi按HB6127—1986第3.9節公式依據動靜壓比值求解(當Ma<1時，PS/PT=(1+0.2Ma2)-3.5)。

在單獨一個測壓點的壓力系數不具備明顯的隨來流角變化的規律的情況下，接著研究不同馬赫數下3個測壓點間壓力系數的差值隨來流角變化的規律，結果如圖6、圖7所示。

圖6 Ma=0.4時各測壓點間壓力系數差

圖7 Ma=0.8時各測壓點間壓力系數差

觀察圖6和圖7內各測壓點間壓力系數差隨來流角變化的曲線，可以發現CpL3-L1在隨來流角變化的過程中呈一條曲線，不具備明顯的規律性；而CpL4-L1和CpL4-L3在對應的角度范圍之內，即CpL4-L1在來流角處于0°～60°之間時以及CpL4-L3在來流角處于-30°～30°之間時，與來流角呈近似線性相關的關系。所以根據上述描述可以按照式(12)、式(13)完成探頭本體氣動建模設計和分析，處理得到左、右探頭兩測壓孔間壓力系數之差與局部來流角θ的關系曲線。

ΔCpL4-Lj(MaL,θL)=CpL4(MaL,θL)-CpLj(MaL,θL)

=(PL4-PLj)/QLC=KθL+θ0

(12)

ΔCpR4-Rj(MaR,θR)=CpR4(MaR,θR)-CpRj(MaR,θR)

=(PR4-PRj)/QRC=KθR+θ0

(13)

式中：當θ≤30°時，j=3；當θ>30°時，j=1；K和θ0按氣動樣本數據處理得到，隨Ma不同取具體值。

但在后續算法建模中使用上述關系式時可以發現，本節建立的半圓柱體探頭本體的氣動模型，是建立探頭上4個測壓點間的壓力系數差與半圓柱探頭所在安裝位置處機身局部來流角的之間的映射關系，即在這個模型中，局部來流角作為關系式中的因變量不能作為上式中j取值的判斷依據，應采取Cp1、Cp2、Cp3、Cp4或者CpL4-L3和CpL4-L1作為j取值的判斷依據，由于不同馬赫數下，測壓點間壓力系數差變化規律極為類似，故以Ma=0.4時的情況為例進行分析，仔細分析圖6可以看出，在CpL4-L1≥-1時，CpL4-L1與局部來流角θL呈線性關系；而在CpL4-L1<-1時，CpL4-L3與局部來流角θL呈線性關系。故綜合考慮之后，將CpL4-L1≥-1作為式(12)和式(13)中j值的取值判斷依據，即在式(12)和式(13)中，當CpL4-L1≥-1時，j=1；當CpL4-L1<-1時，j=3；K和θ0按氣動樣本數據處理得到，隨Ma不同取具體值。

將式(10)帶入到式(12)和式(13)中，則可得到探頭所在處局部來流角指示值的表達式：

θLi=[(PL4-PLj)/QLCi-θ0]/K

(14)

θRi=[(PR4-PRj)/QRCi-θ0]/K

(15)

至此，根據探頭上4個測壓點感受到的壓力值即可求解得到探頭所在位置處的局部動壓指示QLCi、局部靜壓指示值PLSi、局部馬赫數指示值MaLi和局部來流角指示值θLi共4個指示值。接下來再進行迭代求解計算，使各參數收斂至系統要求的精度，即可求得左側探頭所在位置處的局部動壓QLC、局部靜壓PLS、局部馬赫數MaL和局部來流角θL。同理可解得右側半圓柱探頭所在位置處的局部動壓QRC、局部靜壓PRS、局部馬赫數MaR和局部來流角θR。

2.3.2 探頭所在位置局部流場氣動模型

本節要完成半圓柱探頭所在位置局部流場氣動模型設計分析工作。前文已述，按照半圓柱探頭大氣數據系統方案構型，半圓柱探頭大氣數據系統包含對稱安裝于前機身左右表面的4個半圓柱體探頭(兩兩對稱)，但在本節進行半圓柱探頭所在位置局部流場氣動模型時，只選擇一對對稱安裝的半圓柱體探頭作為研究對象即可。

首先研究探頭所在位置局部流場的總、靜壓誤差壓力系數模型，左、右靜壓誤差壓力系數隨迎角變化的曲線如圖8所示。

圖8中，圖8(a)和圖8(b)是馬赫數為0.4的情況下，局部流場左右靜壓誤差壓力系數隨迎角變化的曲線，圖8(c)和圖8(d)則是馬赫數為0.8的情況下的曲線。可以明顯看出，不同馬赫數下，其局部流場左、右靜壓誤差壓力系數的變化趨勢幾乎相同，所以在研究其變化規律時，只以馬赫數為0.4的情況為例進行分析。而對比圖8(a)和圖8(b)兩圖可以發現，在側滑角為5°的時候，局部流場左靜壓誤差壓力系數的變化曲線和側滑角為-5°時的局部流場右靜壓誤差壓力系數變化曲線幾乎重合；類似地，側滑角為-5°的時候，局部流場左靜壓誤差壓力系數的變化曲線和側滑角為5°時的局部流場右靜壓誤差壓力系數變化曲線也幾乎重合；而側滑角為0°時的變化曲線則處于側滑角為5°和-5°時變化曲線中間，此時考慮取局部流場左、右靜壓誤差壓力系數的均值，做出圖像如圖9所示。

圖8 不同馬赫數下左右靜壓誤差壓力系數隨迎角變化曲線

圖9 局部流場左右靜壓誤差壓力系數均值隨迎角變化曲線

觀察圖9可以發現，此時3條曲線幾乎重合，而且都貼近側滑角為0°時的曲線，說明局部流場左、右靜壓誤差壓力系數的均值與側滑角的變化沒有較大的關系，只與馬赫數和迎角變化有關。類似地，局部流場左、右總壓誤差壓力系數的變化規律也有一樣的特性，如圖10所示。

圖10 不同馬赫數下局部流場左右總壓誤差壓力系數隨迎角變化曲線

圖11為左右總壓誤差壓力系數均值隨迎角變化曲線。從圖11中可以看出，局部流場左右總壓誤差壓力系數均值的變化規律與局部流場左右靜壓誤差壓力系數均值的變化規律類似，局部流場總壓誤差壓力系數的變化與馬赫數和迎角變化有關，同時與側滑角的絕對值有關。以探頭所在位置的局部流場總、靜壓誤差壓力系數樣本數據集為輸入,可以按照以下定義和式(16)、式(17)完成探頭所在位置局部流場總、靜壓氣動建模設計和分析。

圖11 不同馬赫數下左右總壓誤差壓力系數均值隨迎角變化曲線

(1)取左、右探頭局部流場靜壓均值為指示靜壓PSi，則靜壓誤差壓力系數為

CpS(Ma,αt,βt)=(PSi-PS)/QC

=[CpLS(Ma,αt,βt)+CpRS(Ma,αt,βt)]/2

(16)

(2)取左、右探頭局部流場動總均值為指示總壓Pti，則總壓誤差壓力系數為

Cpt(Ma,αt,βt)=(Pti-Pt)/QC

=[CpLt(Ma,αt,βt)+CpRt(Ma,αt,βt)]/2

(17)

確定總壓誤差壓力系數的變化規律和靜壓誤差壓力系數的變化規律之后，要進一步分析左、右局部來流角θL、θR和迎角、側滑角及馬赫數之間的映射關系。

通過對CFD流場數據進行分析可以發現，在同一迎角、側滑角狀態下，馬赫數改變但θL和θR值保持不變，所以可以認為θL和θR均只為αt、βt兩個變量的函數。此時以Ma=0.4的情況為例，做出在不同側滑角狀態下，θL和θR隨迎角變化的曲線，如圖12、圖13所示。

圖12 θL隨迎角變化曲線

圖13 θR隨迎角變化曲線

對比圖12和圖13中的曲線可以發現以下幾點：

① 不同側滑角情況下的局部來流角隨迎角變化的曲線，其斜率都是相同的；

②θL在側滑角為5°的曲線和θR在側滑角為-5°的曲線基本重合，類似地，θL在側滑角為-5°的曲線和θR在側滑角為5°的曲線基本重合；

③ 在同一組曲線中，以θL的變化曲線為例，側滑角分別為5°和-5°時的θL曲線，基本以側滑角為0°時的曲線為對稱軸，呈對稱關系。

基于上述特點，可以發現，若是求得同一側滑角情況下的θL和θR均值，應為3條基本重合的曲線，這說明，迎角的變化可以只與θL和θR的均值有關，與馬赫數和側滑角均無關，曲線如圖14所示。

圖14 不同側滑角情況下θL和θR均值隨迎角變化曲線

綜上所述，根據θL值(θR值)對應αt、βt兩個變量的二維樣本數據表格，可以按照式(18)～式(20)處理得到相應關系曲線。

(3)左、右局部來流角與真迎角αt關系：

θL(αt，βt)=K1αt+αL0

(18)

θR(αt，βt)=K1αt+αR0

(19)

(4)左、右局部來流角均值與真迎角αt關系：

[θL(αt，βt)+θR(αt，βt)]/2=K1αt+α0

式中：K1、αL0和αR0按樣本數據處理得到，其中K1不隨Ma數和βt變化，為定值，αL0和αR0隨βt不同取值不同。

通過上面的分析已經驗證了θL+θR的均值的變化只與迎角的變化有關，與側滑角還有馬赫數無關，接下來進一步分析θL和θR的差θL-R與迎角、側滑角之間的關系。左右局部來流角差值在不同馬赫數下隨側滑角變化的情況如圖15所示。

圖15 不同迎角下θL-R隨側滑角變化曲線

分析圖15可以發現，在不同的馬赫數和迎角狀態下，θL-R與側滑角呈線性關系，說明θL-R與側滑角存在一種斜率與迎角有關的線性關系，據此可以得到式(20)。

(5)左、右局部來流角之差與βt關系:

ΔθL-R(αt，βt)=θL(αt，βt)-θR(αt，βt)=K2βt

(20)

式中：K2按樣本數據處理得到，隨αt和Ma不同取值不同，但是根據圖15可以看出，K2的取值并不是隨αt以及Ma線性變化，所以最終K2的取值可以隨著αt和Ma的變化進行線性插值計算。

綜上所述，可以得到如下公式。

指示動壓：

QCi=(QLC+QRC)/2

(21)

指示靜壓：

PSi=(PLS+PRS)/2

(22)

指示總壓：

Pti=PSi+Qci

(23)

指示馬赫數：

Mai=f(QCi/PSi)

(24)

真迎角：

αt=[(θL+θR)/2-α0]/K1

(25)

真側滑角：

βt=(θL-θR)/K2

(26)

與2.3.1節內容類似，至此按照上述描述，根據左右局部來流角、局部總/靜壓，即可求得遠端自由流場的指示靜壓PSi、指示動壓QCi、指示總壓Pti、指示馬赫數Mai四個指示值和真迎角αt、真側滑角βt。接下來進行迭代求解計算，求得真實靜壓PS、真實動壓QC、真實馬赫數Ma。

2.3.3 逆算法模型

由于仿真計算驗證的需要設計了逆算法，逆算法是基于同樣的氣動建模設計輸入條件，由PS、QC、Ma、αt、βt計算得到左、右探頭8個測壓點壓力。基本思路是通過對Ma、αt、βt進行查表和插值計算，得到左右探頭處局部流場信息，再由局部流場的局部總壓、局部靜壓、局部馬赫數和局部來流角信息進一步通過模型公式和插值計算得到半圓柱探頭上4個測壓點處的壓力。

3 算法仿真驗證與分析

3.1 仿真結果

針對算法仿真驗證的基本思路是，截取一段真實的飛參數據作為輸入，以動壓QC、靜壓PS、馬赫數Ma、迎角αt、側滑角βt為輸入值，通過逆向算法計算得到左右探頭上各4個測壓孔處壓力值。再將8個測壓孔處壓力值帶入到半圓柱探頭大氣數據系統解算模型中來，得到解算出的解算動壓QC′、解算靜壓PS′、解算馬赫數Ma′、解算迎角αt′、解算側滑角βt′，對比解算前后對應數據的差值來驗證算法模型的準確性。

仿真的輸入數據截取了一段某飛機試飛的飛參數據，在該段時間內，飛機呈先上升再平飛，先加速再勻速的飛行狀態，迎角和側滑角均有較大的變化幅度，具體的仿真結果如圖16～圖20所示。

圖16 迎角解算結果

圖17 側滑角解算結果

圖18 馬赫數解算結果

圖19 靜壓解算結果

從圖16～圖20的仿真結果曲線來看，每組仿真曲線的圖(a)為各大氣參數的輸入值和解算值的對比圖，圖(b)為各大氣參數的輸入值和解算值的差值曲線圖，迎角的誤差值基本控制在0.2°以內，側滑角的誤差值基本控制在0.01°以內，馬赫數的誤差值基本控制在0.001以內，靜壓的誤差值控制在30 Pa以內，動壓控制在85 Pa以內。仿真結果表明，本文基于半圓柱體智能探頭設計的大氣數據系統算法模型解算精度達到了預期要求，且迭代收斂穩定，能夠基本滿足工程使用上的要求。

圖20 動壓解算結果

3.2 誤差分析

筆者設計的基于半圓柱體智能探頭的大氣數據系統算法模型，主要誤差來源于在正向建模時，由于樣本數量的限制帶來的數據擬合誤差，在模型設計的后續階段會采取繼續通過CFD計算或風洞試驗增加樣本數量[14]，提高氣動建模精度。

4 結束語

本文設計了一種基于新型半圓柱形智能探頭的大氣數據系統算法模型，介紹了半圓柱探頭外形信息和安裝要求，詳細說明了算法模型的建模思路和過程，并通過仿真驗證證明了該算法模型精度滿足設計要求，迭代收斂穩定。該種模型算法的研究設計方法可以為后續相關設計工作提供參考,具有較高的理論參考價值和實際應用價值[15]。