計算流體力學—周正貴[(P111-P125)+P176]

1、7.3 網(wǎng)格生成圖7.3（a）為直角坐標系下求解域局部網(wǎng)格，按照節(jié)點序號給出網(wǎng)格線對應的貼體坐標數(shù)值，和，。通常可簡單地將節(jié)點序號作為在曲線坐標中的坐標值，比如。這樣在計算域內即得到正方形網(wǎng)格，如圖7.3（b）所示。如果在求解域內網(wǎng)格已分布完成，對應每一節(jié)點的x y坐標值已知，則坐標變換關系可用節(jié)點坐標值的差分表達式表示，比如點雅可比可表示為 圖 7.3（a）求解域內局部網(wǎng)格 圖7.3（b）計算域內局部網(wǎng)格 因此要得到坐標變換關系必須在求解域內分布網(wǎng)格線，即進行網(wǎng)格生成，而后由這些網(wǎng)格線構成曲線貼體坐標。合理的網(wǎng)格線分布不僅對計算精度有直接影響，甚至會影響計算過程的收斂性。因此，網(wǎng)格生成是流場

2、數(shù)值計算中一個比較關鍵的問題。通常要求生成的網(wǎng)格線要平滑，避免有局部太大的扭曲，不同簇網(wǎng)格線要盡可能正交。網(wǎng)格分布要與物理問題本身相匹配，也就是說疏密分布應與物理量變化率相適應。比如在求解黏性繞流時，在壁面附面層區(qū)流動參數(shù)沿壁面法向變化很劇烈，因此沿此方向網(wǎng)格線要密集分布。對于超音速流動問題在激波附近網(wǎng)格線要加密。網(wǎng)格生成主要有三種方法：代數(shù)方法、微分方程方法和保角變換方法。微分方程方法可以處理各種類型的不規(guī)則邊界，具有較強的通用性，因而應用最為廣泛；代數(shù)方法在工程實際中也有一定的應用；保角變換方法由于適用范圍的局限性很大，目前實際應用較少。在此僅介紹代數(shù)方法和微分方程方法。7.3.1 代數(shù)生

3、成方法代數(shù)方法實際上是一種插值方法，下面以葉柵通道內H型網(wǎng)格生成來說明這一方法的基本思想。如圖7.4，如果采用NS方程計算葉柵通道內流場，要求網(wǎng)格在葉片表面附面層加密。圖7.4 葉柵繞流網(wǎng)格在葉片前后緣附近，由于流動參數(shù)變化較劇烈，因而此區(qū)域網(wǎng)格點分布也要適當加密。為此將流向網(wǎng)格線葉片通道內部分、切向網(wǎng)格線歸為一類，此類網(wǎng)格線上節(jié)點分布要求二端密而中間稀。采用下列數(shù)學表達式對此類網(wǎng)格線上節(jié)點進行分布： （7.26） （7.27） （7.28） ， （7.29）分別為線段起始和終止坐標。為此線段上分布的節(jié)點數(shù)。分別為伸展和阻尼因子。當時網(wǎng)格節(jié)點在此線段上均布，當時網(wǎng)格節(jié)點在線段兩端對稱加密，越小

4、，越大，則二端節(jié)點密度越大，如圖7.5所示。圖7.5 網(wǎng)格節(jié)點分示意 將H型網(wǎng)格進口截面至葉柵通道進口和葉柵通道出口至H型網(wǎng)格出口截面流向網(wǎng)格線歸為另一類，此類網(wǎng)格線上節(jié)點分布要求從一端向另一端節(jié)點分布由密而稀。對于此類網(wǎng)格線采用下述方法進行節(jié)點分布（以進口流向網(wǎng)格線為例）： （7.30） （7.31） （7.32） （7.33）分別為中弧線在前緣點坐標及導數(shù)值。為線段上節(jié)點數(shù)，為線段長度，為前緣處相鄰節(jié)點間距離，為第i個節(jié)點坐標值。采用上式分布網(wǎng)格節(jié)點，可做到節(jié)點間距過渡光滑，在前緣點左右網(wǎng)格節(jié)點間距相等。7.3.2 微分方程生成方法微分方程生成方法根據(jù)所采用微分方程類型可分為雙曲型方程、拋

5、物型方程以及橢圓型方程生成方法。其中橢圓型方程生成方法應用較為廣泛，在此作詳細闡述。如果一溫度分布均勻的高溫機翼置于靜止流場中，在穩(wěn)定狀態(tài)下，四周液體由于熱傳導而呈不同的溫度分布，如圖7.6（a）畫出溫度分布等值線，可將這一簇溫度等值線作為網(wǎng)格線（圖中等線）。為求溫度分布首先要給出控制方程。對于傳熱問題，有： 圖7.6（a） 翼型繞流求解域網(wǎng)格示意 圖7.6（b） 翼型繞流計算域網(wǎng)格示意當定常時，于是。在直角坐標系下展開，得或者寫成因此可以采用上述橢圓型方程生成網(wǎng)格線。還可從液體運動角度分析。假設液體在由內邊界HABCD和外邊界GFE構成的“C”型流道內流動，液體不可壓無粘無旋。這時有。如果采

6、用流函數(shù)等值線構成網(wǎng)格線，而勢函數(shù)等值線構成網(wǎng)格線，則有。在直角坐標系下展開，得 （7.34）邊界條件設定為 （7.35a） （7.35b）圖7.6(a)和圖7.6(b)反映了求解與計算域網(wǎng)格對應關系。利用邊界條件（7.35），通過數(shù)值求解式（7.34）即可得到坐標變換數(shù)值關系。但采用式（7.34）在x y坐標系下計算，求解域邊界不規(guī)則，因而需要將式（7.34）式進行變換，轉化成在計算域內計算。利用式（7.19），（7.34）可轉化成 （7.36）其中由于，故方程保持橢圓性。如果求解內邊界HAB、CD、GFE、GH、DE4個邊界形狀確定，即每一個邊界節(jié)點(x，y)坐標已知，則在計算域矩形邊界上

7、每一邊界節(jié)點x y值都為已知。采用式（7.36）變換是均勻的或者稱為自然的，也就是說求解域內網(wǎng)格線分布僅與邊界上的節(jié)點分布有關。事實上，如果加上一些輔助項還能使坐標拉伸或壓縮，進而對網(wǎng)格線分布密度進行調整。下面就來說明這一點?？疾旆匠趟慕馐?，這就建立了坐標變換關系。顯然這是一個線性變換，坐標上的等間距節(jié)點分布變換到x坐標上也是等間距的。如果方程改成它的解是。為了方便起見，設端點條件為。則有：在坐標上進行等間距節(jié)點分布，對應x坐標上已不是等節(jié)點間距了。當時，處的節(jié)點加密；當時則變稀疏，如圖7.7。可見，改變c值可改變網(wǎng)格節(jié)點的疏密分布。如果希望在確定的所對應的附近改變節(jié)點分布的疏密程度，可以采

8、用下列方程：圖7.7 網(wǎng)格節(jié)點分布規(guī)律 （7.37）其中： 這時如果式（7.37）中則對應的處附近網(wǎng)格節(jié)點分布變稀疏；如果則變密集。如果有n個節(jié)點附近網(wǎng)格疏密程度要調整，則可采用下列關系式： （7.38）當比較大時，由于指數(shù)函數(shù)衰減很快，各點間影響很小。將上述方法推廣到二維問題中，得到： （7.39）其中： （7.40a） （7.40b）源項實現(xiàn)對線分布控制，其中（7.40a）中第1項實現(xiàn)對線與相鄰網(wǎng)格線間距的整體控制圖7.8（a）；第2項實現(xiàn)對線關于控制，如圖7.8（b）。源項實現(xiàn)對線分布控制與源項實現(xiàn)對線分布控制相類似。圖7.8 網(wǎng)格分布控制實際求解是在計算域內，因而要采用式（7.39）的

9、逆變換方程： （7.41）7.3.3 壁面處網(wǎng)格正交性分析 在任意線坐標系下，二維非定常可壓縮流NS方程可寫成： （7.42）將任一流動參數(shù)f的x，y偏導數(shù)表示成任意曲線坐標關系： （7.43）其中，J為坐標變換雅可比。如圖7.9，考察壁面節(jié)點。采用一階精度向前差分，采用二階精度中心差分，即圖7.9 近壁區(qū)網(wǎng)格對于葉柵通道流，如果采用H型網(wǎng)格計算，等線（葉柵切向）通常與y軸平行，這時?？赏茖С龅碾x散誤差表達式如下：由于，所以有： （7.44a）又有：。這時方程7.44（a）變?yōu)?（7.44b）同理可得： （7.44c）方程7.44（a）、（b）、（c）表明，流動參數(shù)y方向偏導數(shù)離散誤差與網(wǎng)格線

10、夾角無關，而x方向偏導數(shù)離散誤差與成反比。由于壁面附近流動參數(shù)變化劇烈，其二階導數(shù)值相應較大，因此提高網(wǎng)格在壁面處正交性，可減小離散誤差，提高計算精度。下例為平板紊流附層流場計算，考察網(wǎng)格正交性對計算結果影響。平板長度為0.06m，來流。計算網(wǎng)格由平行于平板和斜交于平板兩簇網(wǎng)格線構成，斜交于平板網(wǎng)格線與平板交角分別為，與平板相鄰網(wǎng)格線距平板約為2.0m。取距離平板前緣處垂直截面上計算結果進行比較。圖7.10（a）、（b）分別為附面層內速度整體與近壁區(qū)局部分布比較。表1為壁面剪切應力，其中經(jīng)驗值為采用經(jīng)驗公式計算值，其中DE表示與經(jīng)驗值相對差。由圖7.10（a）、（b）和表1知，當時，附面層內速

11、度分布及壁面剪切應力與經(jīng)驗公式符合很好。比較圖7.10（a）和（b），角對近壁處速度分布影響更大。近壁處速度分布決定壁面剪切應力，由圖7.10（b）和表1，要準確預測壁面阻力，不應小于。 圖7.10（a）附面層內速度分布 圖7.10（b）非?？拷诿嫣幩俣确植急? 壁面剪切應力實驗9080706050403048.9249.4449.5251.6058.7666.5678.2095.61DE(%)0.01.11.25.520.136.059.595.47.3.4 自適應網(wǎng)格簡介在進行網(wǎng)格生成時要求在流動參數(shù)變化較劇烈的區(qū)域分布密集，比如激波附近、尾流區(qū)、附面層區(qū)等，因此在生成網(wǎng)格前要先估計在流

12、場中哪些區(qū)域流動參數(shù)變化較快。對于黏性繞流可以認定在物面附近的附面層區(qū)流動參數(shù)沿物面法向變化較快；但對于有些流動問題，比如超音速流，流場中有無激波、激波位置如何 往往很難預估。所以，在計算過程中，網(wǎng)格的分布最好能根據(jù)計算出的流動參數(shù)空間變化情況不斷調整，這樣就產(chǎn)生了自適應網(wǎng)格（Adaptive Meches）。調整的方法一般有兩種：根據(jù)計算所得流動參數(shù)值，確定新的網(wǎng)格分布，用插值方法計算出新網(wǎng)格節(jié)點上的參數(shù)值；用動網(wǎng)格方法，網(wǎng)格坐標與流動控制方程聯(lián)立求解。下面就一維問題闡述自適應網(wǎng)格的基本思想。引入權函數(shù)，使 （7.45）上式離散化的形式為 （7.46）其中?？梢娭等〈?，網(wǎng)格密度增加，變小；值

13、取小，網(wǎng)格變稀。引入計算域坐標，使，為常數(shù)，通常取網(wǎng)格點數(shù)，這時等于1。于是，式（7.46）變?yōu)?（7.47）上式中常數(shù)若為，x區(qū)間為(0，L)，對應的區(qū)間為(1，N)。則對應的方程為 對上式積分得 于是有： 因此， （7.48）以上積分中以為自變量，如果以x為自變量也可得類似的積分關系。由于 積分得 和 于是有： （7.49） （7.50a） （7.50b）以上給出的式（7.48）和式（7.50）即可用于生成一維自適應網(wǎng)格。它們之間不同的是式（7.48）中權函數(shù)；而在式（7.50）中。網(wǎng)格分布和流動控制方程求解可以分別進行，即每求若干步控制方程的解（這里假設問題是非定常的或迭代求定常解的問題

14、）之后，再根據(jù)生成網(wǎng)格的方程生成新網(wǎng)格。現(xiàn)在的問題是如何取值。最簡單的方法是取 （7.51）為速度的一階導數(shù)，于是根據(jù)式（7.47），有：這樣在速度梯度較大的區(qū)域網(wǎng)格會自動加密。上式進一步可變成 （7.52）這就是說每一網(wǎng)格間距對應于相同的速度增量，即在坐標上是按等取網(wǎng)格點，如圖7.11（a）所示。這種方法的=缺點是當時（在流動均勻的區(qū)域），網(wǎng)格間距。為此可將權重因子作一改變，取相等相等 （7.53）相等（a） （b） （c）圖7.11 權函數(shù)與網(wǎng)格節(jié)點分布這時有 即 為微段弧長，于是式（7.47）變?yōu)?或 （7.54）也就是說網(wǎng)格間距由曲線等弧長條件確定，如圖7.11（b）所示。權函數(shù)更一般

15、的取法是 （7.55）于是在關系曲線上：故 代入式（7.47），得或 （7.56）即在坐標上按等弧長取網(wǎng)格，如圖7.11（c）所示。當時原來斜率大的地方變得更陡，網(wǎng)格加密程度更大。以上作法也有一些缺點，在有些問題中處曲率很大（即很大），如果希望在這些區(qū)域網(wǎng)格加密度，權函數(shù)可用下式： （7.57）其中k為曲線的曲率，為兼顧斜率和曲率的影響，權函數(shù)可以用 （7.58）或 （7.59）以上是一維自適應網(wǎng)格生成的基本方法，這種方法同樣可推廣到二維和三維自適應網(wǎng)格生成。7.3.5 計算網(wǎng)格生成實例卡門翼型繞流計算網(wǎng)格 采用橢圓型微分方程生成求解域中的網(wǎng)格，要將在計算域內的橢圓型方程（7.41）差分離散。如圖7.12，令即（坐標值采用網(wǎng)格節(jié)點序號），得圖7.12 局部網(wǎng)格示意圖 （7.60a） 7.60b）其中： 方程（7.60a、b）可采用點迭代和線迭代方法進行求解，在此采用黎曼點迭代，得 （7.61a） （7.61b）以下采用上述差分方程生成網(wǎng)格，其求解域流場外邊界為一圓，內邊界為卡門翼型，如圖7.13（a），其型面坐標數(shù)據(jù)見附錄V卡門翼型x y坐標值。 圖7.13（a）初始網(wǎng)格 圖7.13（b） 生成的網(wǎng)格圖7.13（c） 計算域網(wǎng)格 采用數(shù)值方法生成網(wǎng)格首先要給定邊界節(jié)點x y坐標，并且還要給出假定初始時刻內部節(jié)點的x