1、昌git空*學課程設計題目：二維超聲速普朗特-邁耶系數波流場的數值解學院：飛行器工程學院專業名稱：飛行器設計與工程班級學號：07034211學生姓名：李桂平指導教師：文U勇二O一O年十一月第一部分1物理問題簡介：普朗特一一邁耶稀疏波的解析解圖-1中，超聲速流圍繞著一個尖的擴張角膨脹，無數個無限弱的馬赫波組成了稀疏波，在尖角處展開成扇形。扇形稀疏波的波頭與來流方向的夾角1,而2是其波尾與下游方向的夾角。和*2稱為馬赫角，定義為：1工1和-Mai和Ma2分別為上下游的馬赫數。通過稀疏波的流動是等熵流動。當流體通過稀疏波后，馬赫數增加，壓力、溫度和密度降低；圖-1中標明了這些變化趨勢。在中心稀疏波前

2、的流動是均勻的，馬赫數為Ma1，而且流動平行與波前的壁面。稀疏波后的流動是均勻的，馬赫數為Ma2，并且平行于下游的壁面。在稀疏波內，流動參數光滑變化，流線彎曲，如圖-1所示。稀疏波內的流動是二維的，唯一的例外是折角的定點，它是一個奇點，壁面流線的方向在此處有一個突然的變化，而且此處的流動參數也是不連續的。這個奇點對流動的數值解會產生影響。給定超聲速來流條件和拐角處的偏轉角一，下游參數是唯一確定的。對于完全氣體，在稀疏波后的流動有精確的解析解，下面給出這個解。流過中心稀疏波的流動，其解析解取決于簡單的關系式'-'-1式中，f是普朗特一一邁耶函數；-是流動偏轉角。對于完全氣體，普朗

3、特一一邁耶函數是Ma和Y勺函數，定義為一arctanvMa21解析解中如下依次得到。對給定的Ma1，從式(2)計算函數打。然后，對給定的偏轉角9,從式(1)得到f2。用這樣得到f2的值，通過求解式求出Ma2。式是關于Ma2的隱式關系式，需要用試湊法求解波后的壓力溫度和密度都可以由等熵流動關系式:珂(曠1)和狀態方程:得到。借助是式(1)式(5)，中心稀疏波后的流動就完全確定了2問題的提法考慮圖-2所示的物理平面。來流馬赫數為2,來流的壓力、密度、溫度分別為：1.01x105N/m2、1.23kg/m3、286.1K。超聲速流動的擴張角0=5.325。，計算區此時，h=h(x)為此時，h=h(x

4、)為域為：x=0到x=60m，壁面到y=40m，如圖-4所示。計算區域內的壓力、密度、溫度、馬赫數等。擴張角定點位置是x=10m40m40+（X-10）tan040m40+（X-10）tan0（0<x<10m）（10<x<65m）初值線。初值線在x=0處，在位于這條鉛垂線的網格點上，初值Mai=2由來流給定。計算從這APt=l.CHxlO5N/m-Ti=2B6.1KWm條線開始并以？X為步長向下游推進。為更好的解決這個丫八10m圖-2問題，下面就對這個問題的解決辦法提出一些理論上的內容，做好準備，以便更明確這個問題的求解第二部分普朗特一一邁耶稀疏波流場的數值解1控制方程

5、定常二維流強守恒形式的控制方程組可以表示如下的通用形式:8F_tK3xSyF和G為列向量，其中:PUpi?+ppuvpu(e+y)+pupupuvpv2+ppv(e+y)+pv考慮沒有體積力的等熵流動。(6)式中的源項J等于零。把(7)式中的列向量每一個分量記作:F】=pnFa=pu2+pF3=Pv7a7b7cF4=pu|e+號)+pu7d對于完全氣體RT1pe=rT=fy-1y-1p因此可消去式(7d)的e,最后得到yu3+v2F產門Pu+PU7e同樣的，可以得到：Gj=puG=puvG3=pv2+p(V2G4=pvle+1+pv=pv+pv8a8b»Hl8ua+va沿流向推進方法

6、的基本思路。在方程(6)中，將x寫在了導數的左邊，y的導數寫在了右邊?？聪聢D一3,如果沿著位于Xo出的處之線給定流場變量是y的函數，那么沿著這條線可以野求出方程中的G的y方向導數，進豐山f而得到F的x方向導數。再由這些方向導數，就可以得到位于x°+?x處的下一條鉛垂線上的流場變量。按這種方式，可以從沿著初值線給定的流場開始，通過沿x方向以?x為步長的推進得到全流場的解，如圖-2。關于強守恒的試驗表明，數值求解這種形式的方程，會出現一些額外的問題，即：3需要將通量Fi、F2、Fa、F4分解，才能求的原變量；O2向量G的元素G、G、G、G只能用Fi、F2、Fa、F4來表示，而不是像式(8

7、a8d)那樣，用原始變量來表示。下面討論下這兩個問題。對于第一個問題，從通量變量中求出原始變量，結果如下：-BblB-4AC3其中:101112=Fz-Fi以及狀態方程如果用強守恒形式的控制方程進行數值求解，也就是求解方程(6)，直接得到的是通量Fi、F2、F3、F4的值，不是原始變量的值。p、uv、p和T的值必須由式(9)式(13)得出。對于第二個問題，即如何計算方程中的G在給定的網格點上，Fi、F2、F3、F4的值可以直接從方程(6)的數值解得到，所以對于下一個網格點處的計算，將這些值用在G、G、G、G的計算中是有道理的。也就是說，應該用Fi、F2、F3、F4的值直接計算G、G、G、G。而

8、不是先從式(9)式(13)中求的原變量，然后再用這些變量根據式(8ad)計算G、G、G、G。因為，G顯然是F的函數。下面給出這種函數關系。首先，由式(8a)和式(11)，得到-14由方程(7c)和式(8b)，可以得到G，即：G2=F315由式(8c)和式(11)，有-'-116利用式(7b)和式(10)，可以從式(16)中消去p,因為：I-仃將式(17)代入式(16)，得伺+H18最后,G4的表達式如下：19G*=吉PV+卩卑*(兔-班+會(7)2+©"J網格生成與坐標變換為了給流過擴張角的流動建立有限差分解法，必須使用貼體網格系統，如圖-4所示，用xy笛卡爾坐標系

9、統表示的物理平面如圖-4a所示。帶有擴張角的壁面構成了物理平面的下邊界。入流邊界位于x=0，出流邊界位于x=L。上邊界為界。n和E為常數的線組成了計算平面內的有規則的正交網格。在En平面的正交網格上實施有限差分計算，控制流動的偏微分方程組是在轉換后的平面內進行數值求解的。因此，為了能在計算平面內使用他們，必須做適當的變換。即，方程（6）必須變換成E和n表示的形式。導數的變換：上兩式中的度量可以通過變換關系式（20）和式（21）得到，即:dxhaxhth25-4a，用x=E表示擴張角27=3yh式（24）中的一可以可以寫成更簡單的形式，由圖3x的位置，則對x_E:ys=0h=常數對x_E:ys=

10、-(x-E)tan0h=H+(x-E)tan0對這些表達式進行微分，得到：對xE:二一二一故血對x_E:.'dxthh一因此，度量可以表示為：31色（當ME）.26a辦（1-n）（當xAE）.26b將式（22）、式（23）、式（24）、式（25）和式（26）代入式（a）、式（b），就得到了導數的變換：28式（27）中，一由式（26a）或式（26b）給出再來看看方程（6）給出的物理平面內的守恒形式的流動控制方程。由于J=0,方程寫為：0F3G二Sxdy利用式（27）、式（28）對方程（29）進行變換，得到:或di30其中度量一由式（26a）或式（26b）給出。用向量F和G的分量來寫，方程

11、3x31323334（30）就是下面一組方程：連續性方程:X動量方程:丫動量方程:方程（31）到方程（34）就是要在圖-3b所示的計算平面內求解的流動控制方程。2. 推進步長的計算定常無粘超聲速流的控制方程是雙曲線方程，所以沿流向推進求解才是合適的。對于時間推進，根據CFL條件可以得到可允許的最大步長。指出：從物理概念上講，顯式時間推進可允許的的最大時間步長應該小于或等于，聲波從一個網格點運動到相鄰的網格點所需的時間。這種聲波傳播解釋，使我們能夠直觀的確定定長流動的CFL條件。如圖-5中顯示了x站位上垂直排列的網格點。點1處的小擾動沿著該點處的兩條特征線向外傳播。特征線就是流動的馬赫線，它和流

12、線的夾角就是馬赫角卩。假設點1處流線與x軸的夾角為9,那么左行和右行馬赫波與x軸的夾角分別是9卩和9。圖5只給出了點1處的左行馬赫波。設有一條通過點2的水平線，點1的左行特征線與水平線相交與點a。于是，點a和點2之間的水平距離（?x）1為慮;；總35圖咗計當推曲甚的示意圏平面）圖咗計當推曲甚的示意圏平面）36為保證點2處沿流向推進計算穩定性,步長匕的值不能大于軒呀和購認兩根據點2處的CFL條件，為了穩定性，推進步長？x的值不應該大雨？X1；因而點2和點a之間的距離應小于或至多等于，聲波從點1傳到與點2同樣高的位置所需的距離。對于點3處的右行馬赫波，也有類似的結果，設它與過點2的水平線交于點b。

13、點b和點2之間的水平距離伽囂為者中較小的一個。將這種分析應用于Xo處垂直排列的所有網格點，就能給出Xo處沿流向推進的步長？x，表達式為：37上式中二舊-.|.：.小是x=x0處垂直排列的所有網格點上-的絕對中最大的。由于式（20）和（21）定義的坐標變換給出的E=x，那么圖-4b所示的計算區域中，沿流向推進的步長為?E=?x38式中?x由式（37）確定。聯立式（37）和式（38），并引入柯朗數C,得到?E應滿足的穩定條件Aytan(9+|maxAytan(9+|max39CFL條件要求上式中的C_1第三部分最終結果用空間推進計算得到的結果如下:圖-6列出的是x=8.9375m,11.375m,

14、23.5625m,31.6875m,47.7375m和64.1875m處的x方向速度分量u對縱坐標y的函數圖。速度u/(m/s)圖-7列出的是圖-6速度u對縱坐標y的函數x=8.9375m，11.375m,23.5625m,31.6875m,47.7375m和64.1875m處的壓力64.1875m處的壓力p對縱坐標y的函數圖。圖-7壓力p對縱坐標y的函數圖-8列出的是x=8.9375m,11.375m,23.5625m,31.6875m,47.7375m和64.1875m處的溫度T對縱坐標y的函數圖。290T/kIJ2270260250y/(y-1)mx=E9375m-x=11375m3-X=23.5625mx=316S75mx=47.9375in6x=641S75圖-8溫度T對縱坐標y的函數圖-8溫度T對縱坐標y的函數圖-9列出的是x=8.9375m，11.375m,23.5625m,31.6875m,47.7375m和64.1875m處的馬赫數Ma對縱坐標y的函數圖。圖-9馬赫數Ma對縱坐標y的函數第四部分小結通過這次課程設計，對空間推進的原理有了較深的理解。空間推進用的是定常流的守恒方程組。根據求解區域的形狀，使用貼體坐標系。練習了網格生成的某些方法，并使用變換后的控制方程組。還用到了撲捉波的技術。捕捉波應該使用守恒式的控制方程組，還要添加適當的人工粘性