2013年全國研究生數(shù)學(xué)建模競賽A題

1、2013年（第十屆）全國研究生數(shù)學(xué)建模競賽A題變循環(huán)發(fā)動機(jī)部件法建模及優(yōu)化由飛機(jī)/發(fā)動機(jī)設(shè)計原理可知，對于持續(xù)高馬赫數(shù)飛行任務(wù)，需要高單位推力的渦噴循環(huán)，反之，如果任務(wù)強(qiáng)調(diào)低馬赫數(shù)和長航程，就需要低耗油率的渦扇循環(huán)。雙涵道變循環(huán)發(fā)動機(jī)可以同時具備高速時的大推力與低速時的低油耗。變循環(huán)發(fā)動機(jī)的內(nèi)在性能優(yōu)勢，受到了各航空強(qiáng)國的重視，是目前航空發(fā)動機(jī)的重要研究方向。1變循環(huán)發(fā)動機(jī)的構(gòu)、造及基本原理1.1基本構(gòu)造雙涵道變循環(huán)發(fā)動機(jī)的基本構(gòu)造見圖1、圖2，其主要部件有：進(jìn)氣道、風(fēng)扇、副外涵道、CDFS涵道、核心驅(qū)動風(fēng)扇級（CDFS）、主外涵道、前混合器、高壓壓氣機(jī)、主燃燒室、高壓渦輪、低壓渦輪、后混合器

2、、加力燃燒室、尾噴管。雙涵道模式下，選擇活門和后混合器（后VABI）全部打開；單涵道模式下，選擇活圖2雙涵道變循環(huán)發(fā)動機(jī)結(jié)構(gòu)示意圖圖中數(shù)字序號表示發(fā)動機(jī)各截面參數(shù)的下腳標(biāo)各部件之間的聯(lián)系如圖3所示，變循環(huán)發(fā)動機(jī)為雙轉(zhuǎn)子發(fā)動機(jī)，風(fēng)扇與低壓渦輪相連，CDFS、高壓壓氣機(jī)與高壓渦輪相連，如圖3下方褐色的線所示。藍(lán)色的線表示有部件之間的氣體流動連接（圖3中高壓壓氣機(jī)后不經(jīng)主燃燒室的分流氣流為冷卻氣流，在本題中忽略不計）。油控卜制器CDFFi進(jìn)憶址ELE紳輪心混合器典歯4彳卡*-;加力悔燒寶圖3變循環(huán)發(fā)動機(jī)工作原理圖1.2工作原理變循環(huán)發(fā)動機(jī)有兩種工作模式，分別為渦噴模式和渦扇模式。發(fā)動機(jī)在亞音速巡航的

3、低功率工作狀態(tài)，風(fēng)扇后的模式轉(zhuǎn)換活門因?yàn)楦蓖夂c風(fēng)扇后的壓差打開，使更多空氣進(jìn)入副外涵，同時前混合器面積開大，打開后混合器，增大涵道比，降低油耗，此時為發(fā)動機(jī)的渦扇模式。發(fā)動機(jī)在超音速巡航、加速、爬升狀態(tài)時，前混合器面積關(guān)小，副外涵壓力增大，選擇活門關(guān)閉，迫使絕大部分氣體進(jìn)入核心機(jī)，產(chǎn)生高的推力，此時為發(fā)動機(jī)的渦噴模式。變循環(huán)發(fā)動機(jī)部件建模法燃?xì)鉁u輪發(fā)動機(jī)的特性可以用實(shí)驗(yàn)方法和計算方法獲得。但實(shí)驗(yàn)的方法需要研制復(fù)雜的設(shè)備、投入巨額的資金和消耗巨大的能源，因此實(shí)驗(yàn)的方法不可能經(jīng)常采用。隨著計算能力的不斷提高，發(fā)動機(jī)數(shù)學(xué)模型研究的不斷深入，計算機(jī)仿真精度也在不斷提高，一定程度上彌補(bǔ)了實(shí)驗(yàn)方法的不足

4、，尤其是在發(fā)動機(jī)型號研制過程中，燃?xì)鉁u輪發(fā)動機(jī)計算機(jī)仿真技術(shù)發(fā)揮了不可替代的作用。燃?xì)鉁u輪發(fā)動機(jī)是由進(jìn)氣道、壓氣機(jī)、主燃燒室、渦輪、噴管等部件組成的。如果計算機(jī)能夠?qū)@些部件的性能進(jìn)行準(zhǔn)確的模擬，那么也就能準(zhǔn)確地模擬整個發(fā)動機(jī)的性能。這種建立在準(zhǔn)確模擬發(fā)動機(jī)各部件性能基礎(chǔ)上的發(fā)動機(jī)性能計算方法，稱為部件法。該方法是建立在發(fā)動機(jī)各部件特性已知的基礎(chǔ)上的，因此是計算精度較高的一種方法。附錄1分別對變循環(huán)發(fā)動機(jī)每個部件的計算公式進(jìn)行了逐一介紹。發(fā)動機(jī)平衡方程發(fā)動機(jī)各部件匹配工作時，受如下7個平衡方程制約。1）低壓軸功率平衡NN耳=0CLTLmL其中N是風(fēng)扇消耗功率，N是低壓渦輪發(fā)出功率率2）CLTL

5、（1）耳=0.99為中間軸機(jī)械效mL高壓軸功率平衡（2）N是高壓渦輪發(fā)出功率，THN+N-N耳二0CHCDFSTHmHNCHmH和N分別是高壓壓氣機(jī)和CDFS的消耗功率,CDFS=0.99是高速軸的機(jī)械效率。WW=0g41g413）高壓渦輪進(jìn)口截面流量平衡3）W是高壓渦輪進(jìn)口截面氣體流量,即主燃燒室出口氣體流量和冷卻空氣流量,g41W是通過高壓渦輪特性數(shù)據(jù)線性插值得到的高壓渦輪流量，這里忽略冷卻的空g41氣流量。4）低壓渦輪進(jìn)口截面流量平衡(4)WW=0g45g45W是低壓渦輪進(jìn)口截面流量，由主燃燒室出口氣體流量和冷卻空氣流量計算得g45到，W是通過低壓渦輪特性數(shù)據(jù)線性插值得到的低壓渦輪流量，

6、這里忽略冷卻g45的空氣流量。5）后混合器靜壓平衡pp=0（5）6162p和p分別是后混合器內(nèi)、外涵道(主外涵道)的靜壓，二者應(yīng)平衡。后混合6162器見圖2的6截面.尾噴管面積平衡A-A二0(6)88A為給定的尾噴管8截面的面積，這里A=9.4575e+003、A為按附錄1尾噴管的888有關(guān)公式計算出的尾噴管8截面的面積，二者應(yīng)相等。風(fēng)扇出口流量平衡TOC o 1-5 h zW-W-W二0(7)a2a21a13其中風(fēng)扇出口的流量W分流為副外涵流量W和CDFS進(jìn)口流量W，三者之間a2a13a21應(yīng)存在平衡關(guān)系。副外涵流量W是按附錄1中前混合器的有關(guān)公式計算出的。a13方程(1),(2),(3),

7、(4),(5),(6),(7)中的變量N,N,N,N,N，CLTLCHCDFSTHW,W,W,W,p,p,A,W,W,W各量可分別由附錄1中的公式轉(zhuǎn)化為g41g41g45g4561628a2a21a13表3.1所列12個設(shè)計中需要給出的發(fā)動機(jī)各部件參數(shù)的函數(shù)。表3.1發(fā)動機(jī)參數(shù)說明試給參數(shù)說明試給參數(shù)說明nL低壓轉(zhuǎn)速(風(fēng)扇、低壓渦輪物理轉(zhuǎn)速)ZTH高壓渦輪壓比函數(shù)值nH高壓轉(zhuǎn)速(高壓壓氣機(jī)、CDFS、高壓渦輪物理轉(zhuǎn)速)ZTL低壓渦輪壓比函數(shù)值Z風(fēng)扇壓比函數(shù)值aL風(fēng)扇導(dǎo)葉角CLZCDFS壓比函數(shù)值acdfsCDFS導(dǎo)葉角CDFSZ高壓壓氣機(jī)壓比函數(shù)值aH高壓壓氣機(jī)導(dǎo)葉角CHT*4主燃燒室出口溫度

8、aCH低壓渦輪導(dǎo)葉角要解決的問題請你們完成以下幾個問題：一、1)請畫出附錄4中風(fēng)扇特性數(shù)據(jù)表中流量隨壓比函數(shù)值變化的圖形。2)設(shè)在發(fā)動機(jī)飛行高度H=11km，飛行馬赫數(shù)Ma=0.8的亞音速巡航點(diǎn)，導(dǎo)葉角度均設(shè)置為0。，風(fēng)扇和CDFS的物理轉(zhuǎn)速都為0.95,風(fēng)扇和CDFS的壓比函數(shù)值都為0.5,求風(fēng)扇和CDFS的出口總溫、總壓和流量。二、設(shè)在發(fā)動機(jī)飛行高度H=11km，飛行馬赫數(shù)Ma=0.8的亞音速巡航點(diǎn)，采用雙涵道模式，導(dǎo)葉角度均設(shè)置為0。，選擇活門完全打開，副外涵道面積設(shè)為1.8395e+003，后混合器出口總面積設(shè)置為2.8518e+004，尾噴管喉道面積A二9.5544e+003，n=0

9、.85。請運(yùn)用或設(shè)計適當(dāng)?shù)乃惴ㄇ蠼庥砂l(fā)動機(jī)7個平衡方8L程(1)，(2)，(3)，(4)，(5)，(6)，(7)組成的非線性方程組。要求陳述算法的關(guān)鍵步驟及其解釋，盡可能討論算法的有效性。如果你們隊(duì)還有時間，請研究下面的問題：三、1）設(shè)在發(fā)動機(jī)飛行高度H=11km，飛行馬赫數(shù)Ma=1.5的超音速巡航點(diǎn)，發(fā)動機(jī)采用單涵道模式，將選擇活門面積設(shè)置為0，風(fēng)扇導(dǎo)葉角度、高壓壓氣機(jī)導(dǎo)葉角度、高壓渦輪導(dǎo)葉角度均設(shè)置為0。，后混合器面積設(shè)置為2.8518e+004。請問發(fā)動機(jī)CDFS導(dǎo)葉角度、低壓渦輪導(dǎo)葉角度和噴管喉道面積3個量為多少時，發(fā)動機(jī)的性能最優(yōu)?2）試研究發(fā)動機(jī)飛行高度H=11km，飛行馬赫數(shù)從M

10、a=1.1變化到Ma=1.6,發(fā)動機(jī)特性最優(yōu)時，CDFS導(dǎo)葉角度、低壓渦輪導(dǎo)葉角度，尾噴管喉道面積隨飛行馬赫數(shù)的變化規(guī)律。此時發(fā)動機(jī)采用單涵道模式，將選擇活門面積設(shè)置為0，風(fēng)扇導(dǎo)葉角度、高壓壓氣機(jī)導(dǎo)葉角度、高壓渦輪導(dǎo)葉角度均設(shè)置為0。，后混合器出口總面積設(shè)置為2.8518e+004，后混合器內(nèi)、外涵道面積可調(diào)（即不受附剝后混合器給定的內(nèi)、外涵道面積值的約束）。注：壓比函數(shù)值的定義見附錄3。為了簡單，題中各量的單位不需要轉(zhuǎn)換，直接運(yùn)算認(rèn)為是合理的。附錄1發(fā)動機(jī)部件計算公式附錄2工質(zhì)熱物理性質(zhì)參數(shù)附錄3氣動函數(shù)及其他常用公式附錄4數(shù)據(jù)參考文獻(xiàn)茍學(xué)中，周文祥，黃金泉，變循環(huán)發(fā)動機(jī)部件級建模技術(shù)，航空

11、動力學(xué)報，2013,28（1）：104-111.A題發(fā)動機(jī)評審意見（初稿）一、本題的計分方法，總分110分。1）摘要、寫作10分；2）第一問共30分，其中第1小問5分，第2小問25分；在第2小問中，風(fēng)扇出口的總溫、總壓、流量計算正確且分析敘述清楚18分，CDFS出口的總溫、總壓、流量計算正確且分析敘述清楚7分。3）第二問45分；4）第三問25分，第1小問20分，第2小問5分。第1小問正確給出模型給10分，模型求解10分。二、第一問求解說明本問求解過程中，下面3點(diǎn)在論文計算中應(yīng)有反映。1）本題各部件的模型已給，但各部件之間的聯(lián)系，要靠題意和自己查資料來建模。對于各部件的出口總溫、總壓、流量的分析

12、是貫穿始終的。如文中應(yīng)分析出：風(fēng)扇的出口總溫、總壓等于CDFS的進(jìn)口總溫、總壓。而流量的分析稍微復(fù)雜一點(diǎn)，如風(fēng)扇的出口的流量等于CDFS進(jìn)口流量和副外涵道進(jìn)口流量的和等。2）對特性數(shù)據(jù)的線性插值方法。這里雖然是線性插值，但二維的插值的方法仍然較多，文中應(yīng)敘述插值模型或算法步驟，有分析比較更好如果是用Matlab軟件進(jìn)行插值，本題使用griddata命令可直接完成插值，如果使用interp2和interpl必須要輔助于其他的方法，否則不能完成線性插值。3）利用熵函數(shù)和焓函數(shù)反求溫度?？梢允褂枚址?，也可以使用別的求解方法。如果使用Matlab求解的應(yīng)注明求解命令。如fzero（fun,xO）等。

13、注1：這三點(diǎn)在以后各問中也必須考慮，但是在論文的的敘述過程中不一定反映；如果這三點(diǎn)出現(xiàn)明顯的錯誤則以后各問的求解結(jié)果也一定是錯的。因此上在打總分時要考慮這個細(xì)節(jié)。注2：由于采用的插值方法和求解方程的方法不同，造成了求解的結(jié)果不同，但在本問的求解結(jié)果中差異應(yīng)該不大。三、對第二問求解及評分的說明本問因求解方法靈活，因此給分細(xì)節(jié)暫時不好確定，還請各位專家多討論。1）在本問求解中除了整體上對平衡方程的求解外，還應(yīng)包括對氣動方程的分析和反解，但是研究生們可能在程序中有考慮，但不一定寫進(jìn)了論文。2）燃燒室后各部件的熵、焓函數(shù)應(yīng)使用燃?xì)獾撵?、焓?shù)即要考慮油氣比；氣體常數(shù)，也應(yīng)該同燃燒室前面的系數(shù)不同；前混合

14、器的流量系數(shù)應(yīng)同后混合器的流量系數(shù)不同。3）平衡方程（7）相對難于理解，題中已經(jīng)做了提示，風(fēng)扇、CDFS的流量是用其特性數(shù)據(jù)線性插值出來的，而副外涵道流量是用前混合器計算的前四步求出來的。4）氣動函數(shù)的反解方法可以有多種，如二分法，牛頓迭代法等，也可以直接使用Matalab工具箱求解。5）平衡方程的求解方法，傳統(tǒng)方法使用N+1參量法，其缺點(diǎn)是初始值不好選取，如果選的不合適，就有可能得不到解。在本題我們希望看到更多的求解方法，如粒子群算法等，對初值選取方法的討論是必須的。但不管使用什么解法，都要結(jié)合本問，寫出明確的算法思想或算法步驟等。6）算法的有效性分析等。題目中已經(jīng)提出盡可能的討論算法的有效

15、性，在評分中應(yīng)該有所反應(yīng)。7）平衡方程求解中，應(yīng)使用相對誤差來度量殘差誤差，因?yàn)楦髁恐g的單位不同，相差量級不同。8）計算應(yīng)有結(jié)果，平衡方程應(yīng)敘述殘差的誤差級。本問很難給出一個統(tǒng)一的求解結(jié)果，因?yàn)槌踔挡煌逯捣椒ú煌?，氣動函?shù)、熵函數(shù)、焓函數(shù)的反解方法不同，其結(jié)果就不相同。本問中各平衡方程的殘差的量級不一樣，特別是平衡方程（7）的相對殘差級在0.6左右，其余大部分平衡方程的相對殘差不超過0.001個量級。9）本題的第二問求解誤差比較大，因此很自然的想法是調(diào)整后混合器的面積來提高精度，但是由于平衡方程（7）不涉及后混合器，所以調(diào)整后混合器的面積不能改進(jìn)平衡方程（7）的求解誤差?；谝陨峡紤]第二

16、問的計分方法如下：1）平衡方程求解算法敘述正確10分。2）平衡方程求解算法有效分析10分。3）使用已給的算法對平衡方程的求解分析正確，求解合理（包括各部件的求解的分析，部件之間聯(lián)系的討論、氣動函數(shù)反解的討論、使用相對誤差來度量殘差誤差，算法初值的討論等）25分。四、第二問的一種求解方法1）.第二問的求解除了考慮7個平衡方程。還應(yīng)該考慮q）,（九）小于1這個條件。這樣解7個平衡方程，實(shí)際上還應(yīng)考慮如下的兩個不等式，即WT*g1251251（1）（kp*A）m125125-3351（2）p*225而這兩個不等式只涉及風(fēng)扇、CDFS、高壓氣機(jī)的參數(shù)，低壓轉(zhuǎn)速、高壓轉(zhuǎn)速、和他們的壓比函數(shù)值，在低壓轉(zhuǎn)速

17、給定的條件下，實(shí)際只有四個參數(shù)，這四個參數(shù)都位于0到1之間?？梢圆捎帽闅v的方式找到滿足不等式的解。然后把滿足不等式的解代入到平衡方程（7）中，計算相對殘差，然而滿足上面兩個不等式的殘差很大，不超過0.6的參數(shù)已經(jīng)很少，可將其代入到其余六個平衡方程中計算相對殘差。注：平衡方程為：W-W-W=0a2a21a13以上求解思路不需要考慮初值，研究生們不一定能夠能想到。五、第三問是一個優(yōu)化模型和控制模型。本問中的優(yōu)化模型對發(fā)動機(jī)的特性最優(yōu)的理解不同，建立的模型也會不相同。發(fā)動機(jī)特性最優(yōu)要考慮的特性有：推力最大、耗油率最低、喘振裕度最小、不超溫、不超轉(zhuǎn)等。優(yōu)化模型的約束條件包括：7個平衡方程，導(dǎo)葉角度的范

18、圍，各部件的壓比函數(shù)值的范圍在0到1之間，主燃燒室的出口溫度在0到2000，高、低壓轉(zhuǎn)速位于0到1之間，尾噴管喉道面積大于0。所建立的模型可以是單目標(biāo)優(yōu)化模型，也可以是多目標(biāo)優(yōu)化模型。但發(fā)動機(jī)的最優(yōu)特性至少要考慮推力和耗油兩個特性。由于采用單涵道模式，選擇活門已經(jīng)關(guān)閉，平衡方程(7)自然成立，所以不需要考慮前混合器，平衡方程變的簡單了。本問的第2小問可以在第1小問的基礎(chǔ)上，對模型稍加修改后(如再加上約束條件:飛行Ma數(shù)在1.1到1.6之間，后混合器內(nèi)外涵道的面積和為2.8518e+004且都大于0等)求解，也可以另建模型。六、第一問答案量流110100908070605040300.10.20.30.40.50.60.70.80.91壓比函數(shù)值1)利用附錄4數(shù)據(jù)畫出風(fēng)扇流量隨壓壓比函數(shù)值變化的圖形。200不同的物理轉(zhuǎn)速，流量隨壓比函數(shù)值變化的圖形不同，圖中應(yīng)有9個曲線，它們近似平行。2)在H=11km，Ma=0.8的亞音速巡航點(diǎn)，導(dǎo)流葉片角度均設(shè)置為0，風(fēng)扇的物理轉(zhuǎn)速都為0.95，風(fēng)扇0.5,H=11.0;Ma=0.8;氣體絕熱指數(shù)k=1.4;氣體常數(shù)R=287.0;計算進(jìn)氣道出口總溫總壓，即為風(fēng)扇進(jìn)口總溫，總壓T*二T=244.381210p*=pm=0.344710i(1)風(fēng)扇的壓比函數(shù)值zz=05、物理轉(zhuǎn)速n二0.95、壓氣機(jī)導(dǎo)葉角度=0解答。風(fēng)扇出口總壓：p*=