周 驍,張海波,王繼強,程 登
(1.南京航空航天大學(xué)江蘇省航空動力系統(tǒng)重點實驗室,南京210016;2.中航工業(yè)西安航空計算技術(shù)研究所,西安710065)
航空飛行器和發(fā)動機之間的相互作用對飛行器整體的性能影響至關(guān)重要[1-2]。隨著對現(xiàn)代飛行器性能與質(zhì)量要求地不斷提高,不能忽視二者間的系統(tǒng)耦合特性,而飛行器與發(fā)動機的一體化建模對了解系統(tǒng)耦合特征,提高系統(tǒng)性能起著非常重要的作用。如果仍然沿用航空發(fā)動機和飛行器獨立分析、分離建模的方法,則不能在設(shè)計過程中充分考慮二者的耦合特性,從而使得系統(tǒng)整體性能得不到充分發(fā)揮[3-4]。因此,需要建立發(fā)動機/飛行器一體化模型來模擬其動靜態(tài)特性,并使其具有實時仿真能力,從而為綜合控制系統(tǒng)設(shè)計提供模型支持,使飛/發(fā)一體化盡可能地發(fā)揮其潛在性能[5-6]。
目前,發(fā)動機/飛行器模型大部分是由C語言編寫的部件級模型,該平臺缺少控制和分析方面的工具,因此需要建立1個可以快速實現(xiàn)控制系統(tǒng)分析、診斷的圖形化仿真平臺[7-8]??焖倏刂圃停≧apid ControlPrototype,RCP)技術(shù)是國外近幾年發(fā)展成熟的1種仿真技術(shù),其原理是在控制系統(tǒng)開發(fā)的初期階段,快速地建立控制器或被控對象模型,并通過對整個控制系統(tǒng)進行多次在線或離線試驗以驗證方案的可行性[9]。歐美等發(fā)達國家已將該技術(shù)廣泛應(yīng)用于航空航天、汽車、軍事裝備等領(lǐng)域,并取得一定成果。如AgrestiM[10]、KharyI[11]建立了基于Simulink的發(fā)動機模型;國內(nèi)王建峰[12]、夏飛[13]也開展了相關(guān)研究,但都僅針對于航空發(fā)動機,而沒有考慮飛機或直升機與發(fā)動機之間的動態(tài)耦合,而此耦合對于直/發(fā)一體化建模與綜合控制具有重要影響。國內(nèi)對實物在回路仿真應(yīng)用方面還存在差距,導(dǎo)致對直/發(fā)一體化模型的綜合控制的仿真試驗研究還不夠充分,因此很有必要開展基于Matlab/Simulink的直/發(fā)一體化綜合模型的研究工作。
本文針對直升機與渦軸發(fā)動機之間的耦合問題,采用Matlab和VC++6.0相結(jié)合的方法,建立基于Matlab/Simulink的直/發(fā)一體化綜合模型,并對其進行數(shù)字仿真驗證。
對于復(fù)雜度極高的直升機/渦軸發(fā)動機模型,若直接使用Simulink內(nèi)置模塊搭建,過程將十分復(fù)雜。本文在直/發(fā)一體化C語言模型[14]的基礎(chǔ)上,建立接口函數(shù)并生成可被S函數(shù)調(diào)用的MEX文件,再通過S-Function模塊調(diào)用生成的動態(tài)鏈接庫,減少了建模的工作量、增加了程序代碼的重用度。
假設(shè)把直升機、渦軸發(fā)動機、控制器全部放入同一S函數(shù)中,將不能實現(xiàn)模塊的通用性。為此本文以分別搭建子模塊后組合的方式建立直升機/渦軸發(fā)動機一體化模型,其流程如下:
(1)通過VC++6.0提取一體化模型的子模塊C語言模型;
(2)在一體化模型的子模塊C語言模型中加入Matlab與C語言的接口函數(shù)——Mexfunction函數(shù),并設(shè)定好輸入、輸出參量的個數(shù)與名稱;
(3)通過Matlab將加入接口函數(shù)Mexfunction的直升機發(fā)動機相關(guān)模塊的C語言模型編譯成動態(tài)鏈接庫(即MEX文件);
(4)“修剪”Matlab自帶的S函數(shù)模塊M文件來調(diào)用生成的動態(tài)鏈接庫(即MEX文件),并另存為工作所需文件名;
(5)運用Simulink內(nèi)置S-Function模塊調(diào)用生成的MEX文件;
(6)根據(jù)各子模塊的輸入、輸出信息,建立直升機/渦軸發(fā)動機一體化模型。
作為非定常、非線性動力學(xué)模型,直升機是直升機/渦軸發(fā)動機一體化模型的重要組成部分,要求其適用于各種飛行狀態(tài),具有較高的置信度。因此,直升機模型以功能劃分分別建立主旋翼、機體、尾槳和尾翼模型。
直升機模塊輸入量包括高度H、前飛速度Vx、前飛指令Vxzl、側(cè)飛指令Vyzl、爬升指令Vzzl、偏航角指令Psizl、旋翼轉(zhuǎn)速指令Ω、旋翼角速度指令Ω˙;輸出量包括總距SM、尾槳總距ST、橫向周期變距A1C、縱向周期變距B1S、前飛速度Vx、側(cè)飛速度Vy、爬升速度Vz。
渦軸發(fā)動機由5級軸流式壓氣機、1級離心式壓氣機、2級軸流式燃氣渦輪以及2級軸流式功率渦輪組成,其結(jié)構(gòu)及特征截面如圖1所示。圖中,0為發(fā)動機未受擾動截面,1為進氣道與渦軸發(fā)動機的相交截面,2、3分別為壓氣機進、出口截面,4為主燃燒室出口截面,44為功率渦輪的進口截面,5、8分別為發(fā)動機排氣裝置進、出口截面。
圖1 渦軸發(fā)動機截面
作為直升機/渦軸發(fā)動機一體化模型的動力裝置,發(fā)動機必須與直升機、控制器進行數(shù)據(jù)交換,并且需要輸出可測參數(shù)用于故障診斷。根據(jù)發(fā)動機原理以及實際工作,選擇渦軸發(fā)動機模塊跟外界交換的變量分別為:輸入量包括高度H、前飛速度Vx、總距SM、扭矩QH、燃油量Wfb、導(dǎo)葉角α;輸出量包括功率渦輪轉(zhuǎn)速NP、燃氣渦輪轉(zhuǎn)速NG、輸出功率HPP、扭矩QH。
控制器對于渦軸發(fā)動機至關(guān)重要。運用ALQR控制算法設(shè)計控制器,該算法具有無窮大的幅值裕度和大于60°的相位裕度,能夠消除穩(wěn)態(tài)誤差,實現(xiàn)對指令的漸進跟蹤[15],因而在多變量魯棒控制中廣泛應(yīng)用。其中渦軸發(fā)動機控制器輸入量包括高度H、前飛速度Vx、總距SM、扭矩QH、功率渦輪轉(zhuǎn)速NP、燃氣渦輪轉(zhuǎn)速NG;輸出量包括燃油量Wfb、導(dǎo)葉角α。
在直升機和渦軸發(fā)動機及其控制器子模塊建模后,可組合成為基于Matlab的直升機/渦軸發(fā)動機一體化模型,正確設(shè)定各子模塊的接收傳遞參量是關(guān)鍵環(huán)節(jié),模塊數(shù)據(jù)傳遞如圖2所示。
在Simulink模型中,子模塊只相當于1個普通的模塊,具有特定的輸入、輸出端口。如果不對直/發(fā)一體化模型進行簡化處理,將影響使用和修改,故需利用Simulink平臺自帶的封裝子系統(tǒng)功能對一體化模型進行封裝處理。
Simulink封裝子系統(tǒng)具有以下優(yōu)點:可以自定義封裝子系統(tǒng)的圖標,便于直觀辨認直升機、渦軸發(fā)動機及其控制器;雙擊封裝后的子系統(tǒng),可彈出參數(shù)對話框,便于更改相關(guān)變量與參數(shù),便于使用;封裝之后的模塊可作為普通模塊添加到Simulink模型應(yīng)用中,也可添加到模塊庫中以供調(diào)用,提高程序代碼的重用率,減少重復(fù)編程;有專有的工作區(qū)域,便于后續(xù)對直升機、渦軸發(fā)動機及其控制器的調(diào)試和修改。
運用Simulink平臺的子系統(tǒng)功能對直升機和渦軸發(fā)動機及其控制器的子模塊進行封裝處理,如圖3所示。從圖中可見,基于Matlab/Simulink的直/發(fā)一體化綜合模型大體上分為直升機、渦軸發(fā)動機及其控制器3個模塊,通過數(shù)據(jù)傳遞建立直/發(fā)一體化建型。
圖2 模塊數(shù)據(jù)傳遞
圖3 基于Matlab/Simulink的直/發(fā)一體化模型
鑒于常規(guī)卡爾曼濾波器存在的缺點,運用改進的卡爾曼濾波器來估計發(fā)動機氣路部件蛻化量。在輸入通道中通過增加估計偏差的積分項,使濾波輸入包含偏差的累計激勵,以消除穩(wěn)態(tài)估計誤差[16]。此外在卡爾曼濾波器的輸入通道之前通過對輸出相似變換至地面狀態(tài),便可使用較少的卡爾曼濾波器通過線性差值估計全包線。改進的卡爾曼濾波器如圖4所示。從圖中可見,以Δy+Δy'代替?zhèn)鹘y(tǒng)的Δy 作為輸入,用以求解發(fā)動機部件蛻化量Δη,其中Δy 為發(fā)動機實際輸出與機載模型輸出的偏差量(經(jīng)相似轉(zhuǎn)換),Δy'為的積分項。
圖4 改進的卡爾曼濾波器
卡爾曼濾波器是基于發(fā)動機小偏差線性模型設(shè)計的,因而首先需對發(fā)動機狀態(tài)變量模型進行設(shè)計,該模型數(shù)學(xué)表達式為
此外,除了部件蛻化引起發(fā)動機工作在非額定狀態(tài)下,在實際工作中的各種噪聲污染以及建模假設(shè)條件等都會給模型帶來誤差,因而設(shè)計中還需考慮到系統(tǒng)誤差和測量誤差的影響。可將式(1)擴展為
式中:w、v 分別為系統(tǒng)、測量噪聲。通常認為w、v 為不相關(guān)的服從正態(tài)分布的零均值白噪聲,其協(xié)方差陣分別為Q 和R。按工程經(jīng)驗選取系統(tǒng)噪聲Q=0.0022I,R=0.0022I,測量噪聲R=0.0022I。
式(2)反映了部件蛻化導(dǎo)致發(fā)動機狀態(tài)量和輸出量的變化。若不考慮外物突然撞擊,發(fā)動機部件的性能蛻化是1個緩慢過程,則可設(shè)Δη˙=0。對于額定發(fā)動機,蛻化量Δη=0;而對于長期服役的發(fā)動機,燃燒室、高低壓渦輪均會不同程度產(chǎn)生性能蛻化,以燃燒室、燃氣渦輪、動力渦輪的效率蛻化為例,設(shè)計發(fā)動機性能蛻化估計模塊,即Δη=[Dcomb,Dgas,Dpow]T。
由于Δη 不能直接獲得,因而在式(2)的基礎(chǔ)上將Δη 增廣為狀態(tài)量,通過估計算法來進一步求取,其增廣形式為
采用卡爾曼濾波器作為狀態(tài)觀測器,可得到如下全維觀測器方程
改進卡爾曼濾波器是在輸入端增加發(fā)動機實際輸出與機載模型輸出的偏差量積分項,用以補償模型與真實發(fā)動機輸出之間的偏差,即卡爾曼濾波器的輸入量為Δy+Δy',由此得到如下改進最優(yōu)狀態(tài)估計方程為
式中:頂標·表示估計值;K=P[C M]TR-1,為卡爾曼濾波器增益陣;P 為如下Riccati方程的解
采用Matlab和VC++6.0相結(jié)合的方法,以分別搭建子模塊而后組合的方式建立基于Matlab/Simulink的直/發(fā)一體化綜合模型。首先運用Matlab與C語言的接口函數(shù)Mexfunction將直升機、發(fā)動機、控制器相關(guān)模塊的C語言模型編譯成動態(tài)鏈接庫,運用Matlab/Simulink內(nèi)置S-Function模塊調(diào)用生成的動態(tài)鏈接庫文件,并封裝確認各模塊的輸入、輸出參量信息,最后根據(jù)接收傳遞參量,建立基于Matlab/Simulink的直/發(fā)一體化綜合模型。本文給出了UH60直升機/T700發(fā)動機綜合模型的設(shè)計算例,在C語言平臺下的UH60直升機/T700發(fā)動機一體化綜合模型與基于Matlab/Simulink的綜合模型相關(guān)參量的變化相應(yīng)曲線分別如圖5、6所示。其中“VC”、“Matlab”分別為C語言、Matlab/Simulink的直/發(fā)一體化綜合模型。
為了驗證該綜合模型的精度,給出了從H=0,Vx=0的狀態(tài)下加速到Vx=40m/s的仿真結(jié)果。從圖5中可見,隨著前飛速度Vx的增加,總距SM與尾槳總距ST緩慢減小,之后趨于平穩(wěn),橫、縱向周期變距A1C與B1S變化相對較小。在前飛速度Vx變化過程中,橫、縱向周期變距A1C與B1S有小幅度變化,隨著前飛速度的穩(wěn)定而趨于平穩(wěn)。從仿真曲線可見,基于Matlab/Simulink與C語言的直/發(fā)綜合模型的變化趨勢基本相符,穩(wěn)態(tài)時曲線基本重合。
從圖6中可見,隨著前飛速度Vx的增大,燃氣渦輪轉(zhuǎn)速NG、燃油流量Wfb、輸出功率HPP、扭矩QH的變化趨勢基本相符,隨著其穩(wěn)定而趨于平穩(wěn)。由于控制器的作用,功率渦輪轉(zhuǎn)速NP基本維持在100%的轉(zhuǎn)速。從動態(tài)過程中的仿真曲線可見,基于Matlab/Simulink與C語言的直/發(fā)綜合模型的變化趨勢基本相符。
綜上所述,由于2種平臺的定時精度、截斷誤差和調(diào)用方式等有所不同,導(dǎo)致在不同平臺下的模型動態(tài)運算稍有不同。但在2種仿真平臺下,動態(tài)仿真曲線趨勢相同,且模型運行到穩(wěn)態(tài)時的仿真曲線重合。因此,基于Matlab/Simulink的直/發(fā)一體化綜合模型在一定程度上滿足了控制系統(tǒng)快速原型設(shè)計的要求,也可作為直/發(fā)綜合優(yōu)化控制、發(fā)動機故障診斷、參數(shù)估計的1種新型仿真平臺。
卡爾曼濾波器是發(fā)動機故障診斷的核心,是基于發(fā)動機小偏差線性模型而設(shè)計的,在發(fā)動機模型上分別對控制量和蛻化量作小階躍,根據(jù)采集到的發(fā)動機模型輸出數(shù)據(jù)擬合求取式(4)的適維矩陣并保存以供調(diào)用。改進的卡爾曼濾波器是在其輸入端增加發(fā)動機(真實)輸出與機載發(fā)動機(額定)模型輸出的偏差量積分項,見式(5)。
圖5 直升機相關(guān)參量對比
圖6 發(fā)動機相關(guān)參量對比
本文以UH60直升機/T700發(fā)動機綜合模型下的仿真結(jié)果為例。在H=500,Vx=30m/s工作狀態(tài)下,分別驗證發(fā)動機單一部件蛻化與多部件同時蛻化的估計效果,分別如圖7、8所示。
圖7 發(fā)動機單部件蛻化估計效果
圖8 發(fā)動機多部件蛻化估計效果
從圖7中可見,在t=10s時,分別設(shè)置發(fā)動機部件蛻化因子Dcomb=0,Dgas=2%,Dpow=0。在5s之內(nèi)可以準確估計出燃氣渦輪效率蛻化量,沒有穩(wěn)態(tài)估計誤差,且Dcomb和Dpow則均接近0,達到了良好的估計效果。同時在t=10s時,分別設(shè)置發(fā)動機部件蛻化因子為Dcomb=0,Dgas=0,Dpow=1%在5s之內(nèi)可準確估計出功率渦輪效率蛻化量,并且估計品質(zhì)良好。
從圖8中可見,在t=10s時,分別設(shè)置發(fā)動機部件蛻化因子為Dcomb=1%,Dgas=2%,Dpow=0。與事先設(shè)置的蛻化因子相比,Dcomb、Dgas、Dpow均達到了良好的估計精度。同時在t=10s,分別設(shè)置發(fā)動機部件蛻化因子為Dcomb=1%,Dgas=0,Dpow=2%,估計效果良好。從仿真效果可見,在基于Matlab/Simulink的直/發(fā)一體化綜合模型仿真平臺上,改進的卡爾曼濾波器對發(fā)動機單一部件蛻化與多部件同時蛻化均能實現(xiàn)快速準確地估計,且估計效果良好。
(1)在考慮直升機與渦軸發(fā)動機相互耦合影響前提下,利用動態(tài)鏈接庫技術(shù)實現(xiàn)了C語言與Matlab的混合建模,建立了基于Matlab/Simulink的UH60直升機/T700發(fā)動機一體化綜合模型,動靜態(tài)精度良好;
(2)在基于Matlab/Simulink的直/發(fā)一體化綜合模型的基礎(chǔ)上,開展了基于改進卡爾曼濾波器的發(fā)動機氣路部件故障診斷研究,且發(fā)動機單部件蛻化與多部件同時蛻化估計效果良好。
[1]SUN Jianguo,Vasilyev V,Hyasov B.Advanced multi-variable control system of aeroengines[M].Beijing:Beihang University Press,2005:10.
[2]方振平,陳萬春,張曙光.航空飛行器飛行動力學(xué)[M].北京:北京航空航天大學(xué)出版社,2005:11.FANG Zhenping,CHEN Wanchun,ZHANG Shuguang.Aviation aircraft flight dynamics[M].Beijing:Beihang University Press,2005:11.(in Chinese)
[3]吳森堂,費玉華.飛行控制系統(tǒng)[M].北京:北京航空航天大學(xué)出版社,2005:358-361.WU Sentang,F(xiàn)EI Yuhua.Flight control system [M].Beijing:Beihang University Press,2005:358-361.(in Chinese)
[4]申安玉,申學(xué)仁,李云保.自動飛行控制系統(tǒng)[M].北京:國防工業(yè)出版社,2003:150-161.SHEN Anyu,SHEN Xueren,LI Yunbao.Automatic flight control system[M].Beijing:National Defense Industry Press,2003:150-161.(in Chinese)
[5]Smith B J,Zagranski R D.Next generation control system for helicopter engines[C]//57th Annual Forum Proceedings of AHS,Washington DC:American Helicopter Society International,2001.
[6]Smith B J,Zagranski R D.Closed loop bench testing of the next generation control system for helicopter engines[C]//58th Annual Forum Proceedings of AHS,Montreal:American Helicopter Society International,2002.
[7]Christopher J.Atkinson.Rapid creation and deployment of software interfaces using Matlab Simulink [C]//AIAA Modeling and Simulation Technologies Conference and Exhibit,Colorado:American Institute of Aeronautics and Astronautics,2006.
[8]Quinn R,Sims J.Improved turbine engine performance,responsiveness, and prognostics using model based control in a hard-ware-in-the-loop simulation[C]//43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,Cincinnati:American Institute of Aeronautics and Astronautics,2007.
[9]張?zhí)旌辏趵^業(yè).微型渦輪發(fā)動機快速原型控制系統(tǒng)[J].航空動力學(xué)報,2007,22(2):274-279.ZHANG Tianhong,WANG Jiye.A rapid prototype control system for micro turbine engine [J].Journal of Aerospace Power,2007,22(2):274-279.(in Chinese)
[10]Agresti M,Csmporeale S M.An objectoriented program for the dynamic simulation of gas turbines[R].ASME 2000-GT-42.
[11]Parker K I,GUO TenHeui.Development of a turbofan engine simulation in a graphical simulation environment [R].NASATM-2003-212543.
[12]王建峰.航空發(fā)動機快速控制原型與實時仿真技術(shù)研究 [D].南京:南京航空航天大學(xué),2011.WANG Jianfeng.Research on rapid control prototyping and realtime simulation for aeroengine[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2011.(in Chinese)
[13]夏飛,黃金泉,周文祥.基于Matlab/Simulink的航空發(fā)動機建模與仿真研究[J].航空動力學(xué)報,2007,22(12):2134-2138.XIA Fei,HUANG Jinquan,ZHOU Wenxiang.Modeling of and simulation research on turbofan engine based on Matlab/Simulink[J].Journal of Aerospace Power,2007,22(12):2134-2138.(in Chinese)
[14]張海波,姚文榮,陳國強.渦軸發(fā)動機/直升機綜合控制仿真平臺設(shè)計[J].推進技術(shù),2011,32(3):383-390.ZHANGHaibo,YAOWenrong,CHEN Guoqiang.Design of a numeric simulation platform for integrated turboshaft engine/helicopter control system[J].Journal of Propulsion Technology,2011,32(3):383-390.(in Chinese)
[15]孫立國,孫健國,張海波.基于直升機/發(fā)動機非線性綜合模型仿真的增廣LQR控制器設(shè)計 [J].航空動力學(xué)報,2010,25(2):471-476.SUN Liguo,SUN Jianguo,ZHANG Haibo. Augmented linear quadratic regulator controller design based on nonlinear integrated helicopter/engine model[J].Journal of Aerospace Power,2010,25(2):471-476.(in Chinese)
[16]張海波,陳霆昊,孫健國.一種新的航空發(fā)動機自適應(yīng)模型設(shè)計與仿真[J].推進技術(shù),2011,32(4):557-563.ZHANG Haibo,CHEN Tinghao,SUN Jianguo.Design and simulation of a new novel engine adaptive model[J].Journal of Propulsion Technology,2011,32(4):557-563.(in Chinese)