旋翼流動(dòng)的塊結(jié)構(gòu)化網(wǎng)格自適應(yīng)方法

肖中云，郭永恒，張 露，崔興達(dá)

(中國空氣動(dòng)力研究與發(fā)展中心，綿陽 621000)

0 引 言

旋翼尾渦是流體在旋轉(zhuǎn)槳葉作用下產(chǎn)生的一種特殊的旋渦流動(dòng)現(xiàn)象，是由槳葉梢部脫出的集中渦和后緣脫出的尾跡面組成的流動(dòng)演化系統(tǒng)，其中槳尖渦在尾跡流動(dòng)中占據(jù)主導(dǎo)作用，是旋渦尾跡流場的骨架[1]。由于槳葉旋轉(zhuǎn)、揮舞等非定常運(yùn)動(dòng)特點(diǎn)，槳尖渦的生成及對后續(xù)槳葉的干擾產(chǎn)生嚴(yán)重的非定常氣動(dòng)載荷，并誘發(fā)槳渦干擾（blade-vortex interaction,BVI）噪聲。當(dāng)前，用CFD 方法模擬旋翼尾渦仍是一大難點(diǎn)[2-4]，通常旋翼槳尖渦直徑約為槳尖弦長的10%，而模擬槳尖渦的網(wǎng)格尺度則可能小到槳尖弦長的1%。由于空間網(wǎng)格尺度不足以匹配旋渦尺寸導(dǎo)致旋渦的數(shù)值耗散過大被抹平，影響到槳渦干擾相關(guān)的氣動(dòng)載荷、噪聲預(yù)測的準(zhǔn)確性。

自適應(yīng)網(wǎng)格加密（adaptive mesh refinement, AMR）是解析局部流動(dòng)特征的最有效方法之一，理論上可以根據(jù)流場特征對網(wǎng)格分辨率進(jìn)行動(dòng)態(tài)調(diào)整，是大范圍內(nèi)模擬旋渦空間演化的理想方法。在傳統(tǒng)計(jì)算方法里面，多塊對接結(jié)構(gòu)網(wǎng)格由于有嚴(yán)格的JKL 指標(biāo)關(guān)系，很難實(shí)現(xiàn)網(wǎng)格的局部加密。非結(jié)構(gòu)網(wǎng)格在數(shù)據(jù)結(jié)構(gòu)上比較靈活，但網(wǎng)格加密后單元質(zhì)量下降及網(wǎng)格信息存儲(chǔ)量大是其面臨的難題。笛卡爾網(wǎng)格由于具有空間網(wǎng)格正交和易于進(jìn)行網(wǎng)格自適應(yīng)加密等特點(diǎn)得到了較好發(fā)展[5-6]，被譽(yù)為是CFD 技術(shù)的未來發(fā)展方向之一。笛卡爾網(wǎng)格自適應(yīng)目前有單元自適應(yīng)和塊自適應(yīng)（block structured AMR, SAMR）兩種方法。單元自適應(yīng)根據(jù)流場計(jì)算需要對網(wǎng)格單元進(jìn)行加密；塊自適應(yīng)[7]指網(wǎng)格加密前后都具有塊結(jié)構(gòu)化的特點(diǎn)，待加密區(qū)域不是單個(gè)的網(wǎng)格單元，而是規(guī)則的立方體區(qū)域。目前發(fā)展的自適應(yīng)加密第三方庫中，采用單元加密方法的如libMesh[8]、P4est[9]等，采用塊自適應(yīng)網(wǎng)格的如CHOMBO[10]、PARAMESH[11]和SAMRAI[12]等。這些三方庫將自適應(yīng)網(wǎng)格加密和并行負(fù)載平衡等函數(shù)封裝為獨(dú)立模塊，將其與所解決的具體物理問題和算法隔離開來，在流體力學(xué)、天文學(xué)、宇宙學(xué)、地質(zhì)學(xué)等多個(gè)領(lǐng)域得到了廣泛應(yīng)用。

塊結(jié)構(gòu)笛卡爾網(wǎng)格結(jié)合了笛卡爾網(wǎng)格與結(jié)構(gòu)網(wǎng)格的優(yōu)點(diǎn)。首先利用笛卡爾網(wǎng)格自適應(yīng)能力強(qiáng)的特點(diǎn)，實(shí)現(xiàn)空間任意區(qū)域的網(wǎng)格自適應(yīng)加密；其次網(wǎng)格具有結(jié)構(gòu)化特征，就使得網(wǎng)格單元的相鄰關(guān)系尤其簡單，并且結(jié)構(gòu)網(wǎng)格的一些高階重構(gòu)計(jì)算方法能夠得到繼續(xù)使用，從而獲得較高的流場計(jì)算效率。本文從上述特點(diǎn)出發(fā)，探索了塊結(jié)構(gòu)笛卡爾網(wǎng)格在旋翼計(jì)算中的生成方法，通過結(jié)構(gòu)貼體網(wǎng)格加背景笛卡爾網(wǎng)格的雙網(wǎng)格技術(shù)構(gòu)建計(jì)算域，發(fā)展了背景網(wǎng)格的幾何自適應(yīng)技術(shù)和基于旋翼尾跡的流場自適應(yīng)技術(shù)，采用該方法在UH-60A 旋翼上構(gòu)建了懸停和前飛流場的自適應(yīng)網(wǎng)格，表明當(dāng)前方法的可行性。

1 塊結(jié)構(gòu)笛卡爾網(wǎng)格生成

1.1 笛卡爾網(wǎng)格劃分方法

塊結(jié)構(gòu)笛卡爾網(wǎng)格劃分方法如圖1 所示。首先定義包圍整個(gè)計(jì)算域的初始長方體網(wǎng)格塊，該網(wǎng)格塊內(nèi)網(wǎng)格為三個(gè)方向均勻正交排列的笛卡爾網(wǎng)格，網(wǎng)格間距分別為（ Δx,Δy,Δz），并且約定每個(gè)方向的網(wǎng)格單元個(gè)數(shù)為偶數(shù)。網(wǎng)格生成是網(wǎng)格塊不斷細(xì)分的遞歸過程，當(dāng)判定一個(gè)網(wǎng)格塊需要加密時(shí)，該網(wǎng)格塊在每個(gè)坐標(biāo)方向上一分為二。這樣經(jīng)過一次網(wǎng)格細(xì)分，二維情況下一個(gè)網(wǎng)格塊被分裂為四個(gè)網(wǎng)格塊，三維情況下一個(gè)網(wǎng)格塊被分裂為八個(gè)網(wǎng)格塊。新產(chǎn)生的子網(wǎng)格塊在三個(gè)方向上的網(wǎng)格維數(shù)同父網(wǎng)格保持不變，網(wǎng)格間距縮小為原來的1/2，當(dāng)父網(wǎng)格塊在某個(gè)方向的網(wǎng)格單元數(shù)為偶數(shù)的情況下，子網(wǎng)格塊在這個(gè)方向上正好兩個(gè)單元對應(yīng)父網(wǎng)格塊的一個(gè)單元。圖1 顯示了笛卡爾網(wǎng)格劃分的前三層網(wǎng)格，其中第三層網(wǎng)格只對第二層的一個(gè)網(wǎng)格塊進(jìn)行了加密，這樣第一、二、三層網(wǎng)格的網(wǎng)格塊數(shù)分別為1、8、8 塊，總的網(wǎng)格塊數(shù)為17 塊。

圖1 塊結(jié)構(gòu)網(wǎng)格劃分方法示意圖Fig. 1 Schematic of block structured grid splitting

根據(jù)上述網(wǎng)格劃分特點(diǎn)，三維笛卡爾網(wǎng)格塊采用八叉樹數(shù)據(jù)結(jié)構(gòu)進(jìn)行管理。如圖2 所示，八叉樹的每個(gè)節(jié)點(diǎn)表示一個(gè)正方體的體積元素，頂部節(jié)點(diǎn)稱為根節(jié)點(diǎn)，也是體積最大的節(jié)點(diǎn)。每個(gè)節(jié)點(diǎn)有八個(gè)子節(jié)點(diǎn)，這八個(gè)子節(jié)點(diǎn)所表示的體積元素加在一起就等于父節(jié)點(diǎn)的體積。沒有后代的節(jié)點(diǎn)稱為葉子節(jié)點(diǎn)，八叉樹葉子節(jié)點(diǎn)代表了分辨率最高的情況。例如將空間某個(gè)區(qū)域的分辨率設(shè)成（ Δxt,Δyt,Δzt），那么覆蓋該區(qū)域的網(wǎng)格塊將不斷細(xì)分，直到每個(gè)葉子網(wǎng)格間距將小于目標(biāo)值為止。八叉樹的葉子節(jié)點(diǎn)覆蓋了整個(gè)計(jì)算域，圖2 顯示的葉子節(jié)點(diǎn)個(gè)數(shù)為15。一般情況下流場計(jì)算只需要在葉子節(jié)點(diǎn)上進(jìn)行，當(dāng)采用多重網(wǎng)格等加速收斂計(jì)算方法時(shí)，中間節(jié)點(diǎn)及根節(jié)點(diǎn)可以參與到計(jì)算中。

圖2 八叉樹數(shù)據(jù)結(jié)構(gòu)示意圖Fig. 2 Schematic of octree data structure

采用上述方式生成的塊結(jié)構(gòu)化笛卡爾網(wǎng)格有很多優(yōu)點(diǎn)。首先是網(wǎng)格描述非常簡潔、存儲(chǔ)量極小，初始網(wǎng)格定義了網(wǎng)格塊中心點(diǎn)坐標(biāo)和三個(gè)方向的尺度及網(wǎng)格間距，新生成網(wǎng)格塊層級和中心點(diǎn)坐標(biāo)由上一級網(wǎng)格塊得到，其他信息如網(wǎng)格點(diǎn)坐標(biāo)、單元體積、單元面積等都可以直接用解析式給出。這樣不需要存儲(chǔ)每個(gè)網(wǎng)格點(diǎn)的信息，從而大大減少對內(nèi)存的需求。其次是和結(jié)構(gòu)化網(wǎng)格一樣，塊結(jié)構(gòu)化的笛卡爾網(wǎng)格流場變量按指標(biāo)順序連續(xù)存放，有利于充分利用高速緩存提高存取效率，并且非常容易在三個(gè)方向上構(gòu)造高階插值模板單元，以很小的成本代價(jià)實(shí)現(xiàn)高階格式。再次，塊結(jié)構(gòu)笛卡爾網(wǎng)格由粗網(wǎng)格不斷細(xì)分得到，子網(wǎng)格和父網(wǎng)格之間構(gòu)成多重網(wǎng)格的嵌套關(guān)系，可以根據(jù)此關(guān)系設(shè)計(jì)多重網(wǎng)格方法，提高流場計(jì)算的收斂速度。

1.2 多層空間填充Z 曲線

當(dāng)前自適應(yīng)網(wǎng)格方法可以簡單地類比于給網(wǎng)格增加或減少補(bǔ)丁，加補(bǔ)丁就是給網(wǎng)格加密，揭補(bǔ)丁就是使網(wǎng)格稀疏。空間填充曲線技術(shù)是為了高效存儲(chǔ)和管理空間網(wǎng)格的一種方法，它的目的是將多維空間數(shù)據(jù)轉(zhuǎn)換到一維空間上，并通過轉(zhuǎn)換后的一維空間索引值存儲(chǔ)和查詢多維數(shù)據(jù)。除了有降低維度的特性外，空間填充曲線還具有數(shù)據(jù)聚類特性，其特點(diǎn)是將空間上鄰近的網(wǎng)格單元映射為一維曲線上盡可能接近的點(diǎn)，因此只需要訪問查詢點(diǎn)的鄰近點(diǎn)，就能夠獲得網(wǎng)格單元的最近鄰單元。常用的空間填充曲線有Z 曲線[13]、Hilbert 曲線和Gray 曲線等，其中Hilbert 曲線和Gray 曲線涉及到方向旋轉(zhuǎn)，映射過程比較復(fù)雜，聚類特性更優(yōu)。相比之下Z 填充曲線的映射過程比較簡單，幾何空間的網(wǎng)格相鄰關(guān)系更容易確定而被網(wǎng)格方法所廣泛采用。本文針對自適應(yīng)網(wǎng)格特點(diǎn)將空間填充Z 曲線進(jìn)行了多層表示，即填充曲線包含了葉子節(jié)點(diǎn)、它們的父節(jié)點(diǎn)及根節(jié)點(diǎn)等全部網(wǎng)格，將所有網(wǎng)格層展開后得到網(wǎng)格關(guān)系如圖3 所示。保留非葉子節(jié)點(diǎn)網(wǎng)格的目的是能夠自由地對網(wǎng)格進(jìn)行加密或者稀疏變換，從圖3 可以看到，通過延拓或者插值計(jì)算，LEVEL 1 和LEVEL 2 的流場信息交換是非常方便和直接的。圖3 同時(shí)給出了表示多層網(wǎng)格的空間填充Z 曲線。其編碼規(guī)則如下，編碼首先從根節(jié)點(diǎn)開始，對于任何一個(gè)節(jié)點(diǎn)，首先判斷其是否存在子節(jié)點(diǎn)，如果有則指向子節(jié)點(diǎn)，如果沒有則指向同一級的鄰居節(jié)點(diǎn)，鄰居節(jié)點(diǎn)之間用Z 曲線連接，當(dāng)子節(jié)點(diǎn)為Z 曲線最后一個(gè)節(jié)點(diǎn)時(shí)回到父節(jié)點(diǎn)。注意這里的節(jié)點(diǎn)不重復(fù)計(jì)數(shù)，如圖3 中當(dāng)出現(xiàn)節(jié)點(diǎn)5 指向節(jié)點(diǎn)1 時(shí)，由于節(jié)點(diǎn)1 已經(jīng)存在于序列中而略過，因此5 的下一個(gè)節(jié)點(diǎn)為6。當(dāng)采用上述方式索引以后，可以看到，任一節(jié)點(diǎn)的鄰居節(jié)點(diǎn)可能是兄弟（相同父節(jié)點(diǎn)），也可能是堂兄弟（父節(jié)點(diǎn)兄弟的孩子），但總的來說，仍具有鄰居節(jié)點(diǎn)位于空間填充曲線附近的聚類特性。

圖3 多層網(wǎng)格的Z 填充曲線Fig. 3 Z filling curves of multi-layer grids

在用空間填充曲線進(jìn)行排序以后，對網(wǎng)格進(jìn)行并行分區(qū)、滿足負(fù)載平衡就變得非常簡單，由于每個(gè)網(wǎng)格塊的網(wǎng)格單元數(shù)相同，因此滿足負(fù)載平衡的并行分區(qū)只需要將一維的空間填充曲線進(jìn)行等量剖分就可以了。當(dāng)然也可以考慮到葉子節(jié)點(diǎn)和非葉子節(jié)點(diǎn)之間計(jì)算量的差異，給不同節(jié)點(diǎn)賦予不同的權(quán)重，采用加權(quán)平均的方法進(jìn)行分區(qū)。圖4 顯示了將21 個(gè)網(wǎng)格塊分配到5 個(gè)進(jìn)程的分區(qū)情況。

圖4 多層網(wǎng)格的并行分區(qū)Fig. 4 Parallel partition of multi-layer grids

1.3 網(wǎng)格加密判則與加密方法

采用一分為八的方式加密網(wǎng)格將導(dǎo)致網(wǎng)格量呈幾何級數(shù)增長。為了讓網(wǎng)格總量可控，本文采取的策略一是限制網(wǎng)格最大層數(shù)，一旦網(wǎng)格層級超過最大值則加密被終止，二是建立合適的加密準(zhǔn)則，將網(wǎng)格加密局限在幾何復(fù)雜、流動(dòng)變化劇烈的局部區(qū)域，并且網(wǎng)格隨流動(dòng)特征變化的自適應(yīng)加密。

限制網(wǎng)格最大最小層數(shù)的作法實(shí)際上是限制了網(wǎng)格單元的最大最小尺寸，最大尺寸規(guī)定的是遠(yuǎn)場的網(wǎng)格尺度，與計(jì)算域的大小相關(guān)；最小尺寸規(guī)定的是模擬流動(dòng)特征的網(wǎng)格最小尺度，比如對于旋翼流動(dòng)來說，網(wǎng)格加密是為了實(shí)現(xiàn)槳尖渦的模擬，因此可以根據(jù)槳尖弦長測算出槳尖渦的渦核直徑，并由此估算出能夠較好捕捉旋渦的網(wǎng)格最小尺度。

網(wǎng)格的局部加密包括了基于幾何特征加密和基于流場特征加密兩種類型。旋翼模擬一般將笛卡爾網(wǎng)格作為背景網(wǎng)格，笛卡爾網(wǎng)格與近場網(wǎng)格之間構(gòu)成重疊關(guān)系，此時(shí)近場外邊界網(wǎng)格尺度就作為重疊區(qū)笛卡爾網(wǎng)格加密的目標(biāo)尺度。流場特征加密主要是基于旋翼槳尖渦的識別與局部加密，目前旋渦識別方法發(fā)展有很多種[14-15],常見的Q-判據(jù)定義方法如下：

其中 Ω表示旋轉(zhuǎn)張量，S表示變形率張量。為了防止網(wǎng)格密度變化過渡劇烈，限定相鄰笛卡爾網(wǎng)格塊的網(wǎng)格尺度相差最大一倍。即當(dāng)對某一網(wǎng)格塊進(jìn)行加密時(shí)，還需要判斷該網(wǎng)格塊的相鄰塊是否存在密度差大于一倍的情況，如果存在，則該相鄰塊也被列入到待加密列表中，以此類推，直到所有塊都滿足條件為止。需要注意的是，這里所謂的相鄰塊包括了面相鄰塊、棱線相鄰塊和角點(diǎn)相鄰塊。對于長方體來說，相鄰塊的數(shù)量最多情況下達(dá)26 個(gè)。圖5 以二維為例對網(wǎng)格加密過程進(jìn)行了說明，圖5（a）顯示了待加密網(wǎng)格塊（用十字虛線表示），同時(shí)識別出周圍三個(gè)網(wǎng)格塊存在密度差大于1 倍的情況（桔色顯示）；圖5（b）對（a）中桔色顯示的三個(gè)網(wǎng)格塊進(jìn)行了加密，同時(shí)又識別出周圍兩個(gè)網(wǎng)格塊存在密度差大于1 倍的情況；圖5（c）為最終的網(wǎng)格加密情況，可以看出每組相鄰塊的網(wǎng)格密度最大相差1 倍；圖5（d）為加密后網(wǎng)格的空間Z 曲線填充情況。圖5（d）只對葉子節(jié)點(diǎn)進(jìn)行了空間填充的Z 曲線顯示，顯示多層網(wǎng)格的Z 填充曲線如圖6 所示。從圖中可以看到，當(dāng)前網(wǎng)格一共包含了5 層，除了最后一層全部是葉子節(jié)點(diǎn)外，其他層凡是有向下箭頭線的均為非葉子節(jié)點(diǎn)，即表明當(dāng)前塊有更細(xì)的網(wǎng)格剖分。圖5（d）可以理解為圖6 在平面內(nèi)的投影，并且刪除非葉子節(jié)點(diǎn)后得到的圖形，這些非葉子節(jié)點(diǎn)不直接參與流場計(jì)算，但在方便網(wǎng)格組織與自適應(yīng)加密中起到作用。

圖5 網(wǎng)格加密過程示意圖Fig. 5 Schematic of mesh refinement

圖6 局部加密網(wǎng)格的Z 填充曲線Fig. 6 Z filling curves for local mesh refinement

2 算例與分析

2.1 旋翼近場貼體網(wǎng)格

為驗(yàn)證網(wǎng)格模型，本文選擇生成UH-60A 旋翼的計(jì)算網(wǎng)格，模擬狀態(tài)包括旋翼的懸停狀態(tài)和前飛狀態(tài)。UH-60A 旋翼共包括四片槳葉，旋翼半徑R=8.178 m，槳葉展弦比15.5，槳葉具有非線性扭轉(zhuǎn)、槳尖后掠等現(xiàn)代旋翼特征。槳葉貼體網(wǎng)格采用多塊對接結(jié)構(gòu)網(wǎng)格（見圖7），其中槳葉剖面拓?fù)浣Y(jié)構(gòu)為O 型，槳葉表面網(wǎng)格為四邊形單元，物面邊界層區(qū)域采用大拉伸比網(wǎng)格進(jìn)行刻畫。

從圖7 中可以看到，旋翼網(wǎng)格分為了4 個(gè)相對獨(dú)立的子網(wǎng)格，由于槳葉在一周旋轉(zhuǎn)運(yùn)動(dòng)過程中，即使是剛性假設(shè)槳葉，也需要進(jìn)行旋轉(zhuǎn)、揮舞、變距等運(yùn)動(dòng)，槳葉與槳葉之間的相對位置會(huì)發(fā)生變化，將單片槳葉網(wǎng)格獨(dú)立作為一個(gè)子網(wǎng)格，這樣就可以讓槳葉網(wǎng)格隨槳葉一起運(yùn)動(dòng)，每個(gè)時(shí)間步的子網(wǎng)格變化通過多次旋轉(zhuǎn)變換得到。近壁區(qū)流場計(jì)算由貼體網(wǎng)格負(fù)責(zé)進(jìn)行，背景網(wǎng)格深入到壁面的部分經(jīng)過重疊挖洞以后被當(dāng)作洞內(nèi)點(diǎn)不參與計(jì)算。

圖7 槳葉貼體網(wǎng)格Fig. 7 Body fitted grids of the blade

2.2 背景網(wǎng)格重疊邊界的幾何自適應(yīng)

在生成近場貼體網(wǎng)格以后，背景網(wǎng)格采用笛卡爾網(wǎng)格生成方法，即采用雙網(wǎng)格技術(shù)構(gòu)建流場計(jì)算網(wǎng)格。兩組網(wǎng)格的計(jì)算域劃分如圖8 所示，考慮到當(dāng)前旋翼半徑R= 8.178 m，定義背景網(wǎng)格的外廓尺寸為x,y,z∈[-50 m, 50 m]的立方體區(qū)間，網(wǎng)格加密層級定義為9 級，這樣最小網(wǎng)格單元尺度約為 Δx=0.048 m。圖中立方體代表的是初始笛卡爾網(wǎng)格塊，塊內(nèi)每個(gè)方向的網(wǎng)格單元個(gè)數(shù)為8 個(gè)。

圖8 貼體網(wǎng)格與背景網(wǎng)格Fig. 8 Body fitted and background grids

背景笛卡爾網(wǎng)格具有能夠自動(dòng)生成的特點(diǎn)，從圖8 所示的初始笛卡爾網(wǎng)格塊開始，不斷對網(wǎng)格進(jìn)行細(xì)分，直到核心區(qū)域網(wǎng)格密度滿足模擬要求為止。在這里背景網(wǎng)格與貼體網(wǎng)格之間構(gòu)成重疊關(guān)系，要求重疊區(qū)域兩組網(wǎng)格的網(wǎng)格尺度相當(dāng)，因此該過程也被稱之為網(wǎng)格幾何自適應(yīng)。由于背景網(wǎng)格要滿足與貼體網(wǎng)格重疊插值的條件，這里以貼體網(wǎng)格的外邊界作為特征對象，對背景網(wǎng)格塊依次進(jìn)行判斷，如果特征對象與網(wǎng)格塊相交或者完全落在網(wǎng)格塊的內(nèi)部，則判定該網(wǎng)格塊為待加密網(wǎng)格塊。對初始網(wǎng)格塊來說，定義網(wǎng)格層級LEVEL 0，在經(jīng)過一次加密后生成8 個(gè)子網(wǎng)格塊，定義為LEVEL 1，以此類推，直到網(wǎng)格塊的網(wǎng)格尺度小于特征對象的網(wǎng)格尺度為止。圖9 顯示的是經(jīng)過幾何自適應(yīng)以后得到的笛卡爾背景網(wǎng)格分布，可以看到經(jīng)過不斷細(xì)分以后，網(wǎng)格密度大的區(qū)域向近場集中；此外網(wǎng)格呈塊狀分布，塊內(nèi)部網(wǎng)格均勻分布，相鄰網(wǎng)格塊的網(wǎng)格間距比值最大為2：1。圖10 顯示了背景網(wǎng)格和貼體網(wǎng)格交界區(qū)域的網(wǎng)格分布?？梢钥吹奖尘熬W(wǎng)格在槳葉根部和尖部都進(jìn)行了加密，在子網(wǎng)格外邊界區(qū)域，加密準(zhǔn)則要求背景網(wǎng)格與貼體網(wǎng)格的網(wǎng)格尺度相當(dāng)，滿足重疊插值對網(wǎng)格的要求。

圖10 近物面區(qū)域的網(wǎng)格分布Fig. 10 Grid distribution in the near wall region

2.3 基于旋翼尾跡的流場自適應(yīng)

除了前面的幾何自適應(yīng)以外，發(fā)展自適應(yīng)網(wǎng)格的主要目的是根據(jù)流場特征對網(wǎng)格進(jìn)行局部加密。為驗(yàn)證網(wǎng)格對旋翼流場的自適應(yīng)能力，本文將通過尾跡模型方法得到尾跡分布作為流場特征，對背景網(wǎng)格進(jìn)行局部加密。

2.3.1 旋翼懸停狀態(tài)

懸停模擬狀態(tài)為旋翼轉(zhuǎn)速Ω= 255 r/m，拉力系數(shù)CT=0.0069。旋翼尾跡用Landgrebe 尾跡模型表示，具體公式如下：

其中：rtip、ztip分別表示槳尖渦的無量綱徑向位置和旋翼軸向位置（用旋翼半徑R無量綱化），Nb表示旋翼槳葉片數(shù)， ψw為槳尖渦尾齡角， θtw為槳葉線性扭轉(zhuǎn)角，其他系數(shù)定義如下：

圖11 為Landgrebe 尾跡模型計(jì)算得到的旋翼懸停狀態(tài)下單片槳葉的槳尖和槳根尾跡分布，可以看到槳尖和槳根渦呈螺旋狀向旋翼下方發(fā)展，槳尖渦的下降速度更快，并且在下降的同時(shí)還向內(nèi)收縮，最終穩(wěn)定在約0.7 倍旋翼半徑位置處。

圖11 懸停狀態(tài)的槳尖和槳根渦線Fig. 11 Root and tip vortex lines of the blade in hover

加密準(zhǔn)則是一旦尾跡線穿過網(wǎng)格塊，則定義該網(wǎng)格塊為待加密網(wǎng)格塊，直到網(wǎng)格密度小于目標(biāo)網(wǎng)格密度為止。這里定義目標(biāo)網(wǎng)格密度為槳尖弦長的2%。圖12 為懸停狀態(tài)下的自適應(yīng)網(wǎng)格分布，可以看到在分布有尾跡線的地方，網(wǎng)格均進(jìn)行了局部加密。由于模擬旋渦的網(wǎng)格尺度要小于槳葉子網(wǎng)格外邊界的平均網(wǎng)格尺度，因此，當(dāng)?shù)氐芽柧W(wǎng)格的分布更密更集中，遠(yuǎn)小于重疊區(qū)附近的網(wǎng)格尺度。

圖12 懸停狀態(tài)的自適應(yīng)網(wǎng)格加密Fig. 12 Adaptive mesh refinement of the rotor in hover

2.3.2 旋翼前飛狀態(tài)

前飛模擬狀態(tài)為UH-60A 旋翼的C8534 前飛狀態(tài)，該狀態(tài)對應(yīng)UH-60A 旋翼的大速度中等過載飛行狀態(tài)：其中拉力系數(shù)CT=0.0069， 槳尖馬赫數(shù)Mtip=0.642 ，前進(jìn)比 μ=0.368， 槳盤傾角αs=-7.31°。

旋翼前飛狀態(tài)尾跡用Beddoes 預(yù)定尾跡模型[16]表示，該模型給出了尾跡分布的代數(shù)表達(dá)式，為研究旋翼尾跡分布以及槳渦干擾BVI 現(xiàn)象提供了方便。Beddoes 預(yù)定尾跡模型中包含了尾齡角 ψw和槳葉方位角 ψb， μ為前進(jìn)比。渦元的軸向位移z則由旋翼軸向速度和局部誘導(dǎo)下洗速度積分組成，則槳尖渦尾跡幾何形狀可表示為：

其中，xtip、ytip、ztip為槳尖渦尾跡的笛卡爾坐標(biāo)，用旋翼半徑R無量綱化，坐標(biāo)原點(diǎn)位于旋翼中心。rv為控制槳尖渦卷起及隨渦齡角向內(nèi)收縮的徑向位置參數(shù)， μ為前進(jìn)比，α為來流迎角， λi為當(dāng)?shù)卣T導(dǎo)入流比，計(jì)算公式如下：

其中： λ0為 平均誘導(dǎo)入流比，E為尾跡傾斜角度參數(shù)。

圖13 給出的是UH-60A 旋翼前飛狀態(tài)下的槳尖尾跡分布三維視圖，由Beddoes 預(yù)定尾跡模型計(jì)算得到。可以看到在前飛狀態(tài)下，旋翼尾跡呈螺旋狀朝旋翼后下方發(fā)展。研究表明，旋翼槳尖渦的渦核直徑為槳尖弦長的1/10 量級，模擬槳尖渦所需網(wǎng)格的長度尺度約為槳尖弦長的1/100 量級，并且槳尖渦的分布范圍很廣，給傳統(tǒng)數(shù)值模擬方法的網(wǎng)格生成帶來了極大挑戰(zhàn)——采用近場全局均勻加密方法導(dǎo)致網(wǎng)格量不可承受，必須發(fā)展網(wǎng)格當(dāng)?shù)丶用艿淖赃m應(yīng)生成技術(shù)。

圖13 前飛狀態(tài)的槳尖渦線Fig. 13 Tip vortex lines of the rotor in forward flight

網(wǎng)格最大加密層數(shù)設(shè)置為9 層，由于網(wǎng)格加密的密度要求高，尾跡線分布范圍廣，最終網(wǎng)格塊數(shù)達(dá)到了4.4×104塊，網(wǎng)格單元數(shù)達(dá)到2.25×107。加密運(yùn)算以單個(gè)進(jìn)程串行方式完成，在主頻3.2G 的Intel I7-8700 CPU 芯片運(yùn)算網(wǎng)格加密時(shí)間約為76 s。圖14是旋翼前飛狀態(tài)下的自適應(yīng)網(wǎng)格分布，可以看到當(dāng)前網(wǎng)格加密方法適應(yīng)了尾跡線的分布情況，在旋翼的后方和下方區(qū)域?qū)W(wǎng)格進(jìn)行了加密，在槳尖尾跡線穿過的地方網(wǎng)格密度達(dá)到最大值。

圖14 前飛狀態(tài)的自適應(yīng)網(wǎng)格加密Fig. 14 Adaptive mesh refinement of the rotor in forward flight

圖15 用分層方式顯示了旋翼附近的網(wǎng)格，位于中間的第6～7 層網(wǎng)格是在前5 層網(wǎng)格上的疊加，位于右邊的第8～9 層網(wǎng)格又是在前7 層網(wǎng)格基礎(chǔ)上的疊加；最大網(wǎng)格密度位于槳尖渦尾跡附近，為尾跡流動(dòng)的精細(xì)模擬創(chuàng)造了條件。

圖15 自適應(yīng)網(wǎng)格的多層結(jié)構(gòu)Fig. 15 Multi-layer structure of adaptive mesh refinement

圖16 給出了旋翼前飛狀態(tài)自適應(yīng)網(wǎng)格的第8 層網(wǎng)格分布，網(wǎng)格最大層級為9 層，第8、第9 層代表了網(wǎng)格的核心加密區(qū)，可以看到這些密網(wǎng)格分布在旋翼附近及旋翼后下方區(qū)域。圖中不同顏色代表網(wǎng)格塊所在的不同分區(qū)，分區(qū)方法利用空間填充Z 曲線的聚類特性，使分在同一個(gè)進(jìn)程的網(wǎng)格塊保持相鄰。同時(shí)需要注意的是由于空間填充Z 曲線在結(jié)束一個(gè)Z 遍歷以后有跳躍的特點(diǎn)，前一個(gè)Z 的最后一個(gè)網(wǎng)格塊與后一個(gè)Z 的第一個(gè)網(wǎng)格塊并不相鄰（見圖5（d）），因此可能出現(xiàn)分在同一個(gè)進(jìn)程的網(wǎng)格塊不相鄰的情況，有研究[17]表明在同一個(gè)進(jìn)程內(nèi)，不相鄰的網(wǎng)格塊最多為兩組，這樣保證了網(wǎng)格塊不過于分散，降低了并行發(fā)送與接收方面的需求。

圖16 塊結(jié)構(gòu)笛卡爾網(wǎng)格的并行分區(qū)Fig. 16 Parallel partition of the block structured Cartesian grids

3 結(jié)束語

本文探討了旋翼計(jì)算中背景笛卡爾網(wǎng)格的生成方法，發(fā)展了塊結(jié)構(gòu)化笛卡爾網(wǎng)格的生成技術(shù)。將該技術(shù)應(yīng)用在UH-60A 旋翼上的懸停和前飛流場，獲得了網(wǎng)格密度分布合理的自適應(yīng)網(wǎng)格，主要結(jié)論如下：

1）采用雙網(wǎng)格建模思路，背景笛卡爾網(wǎng)格由于不需要適應(yīng)物面邊界，網(wǎng)格劃分過程十分快捷高效。背景網(wǎng)格可以通過幾何特征或流場特征對網(wǎng)格進(jìn)行自適應(yīng)加密，使當(dāng)?shù)鼐W(wǎng)格密度滿足重疊插值或流場計(jì)算的要求。

2）塊結(jié)構(gòu)自適應(yīng)笛卡爾網(wǎng)格采用多層網(wǎng)格結(jié)構(gòu)，高效實(shí)現(xiàn)了網(wǎng)格的加密和稀疏操作。用八叉數(shù)數(shù)據(jù)結(jié)構(gòu)行管理網(wǎng)格塊，相對于傳統(tǒng)笛卡爾網(wǎng)格管理網(wǎng)格單元而言，整個(gè)八叉數(shù)所管理的節(jié)點(diǎn)數(shù)大大降低，提高了整個(gè)數(shù)據(jù)結(jié)構(gòu)的管理效率。

3）采用多層空間填充Z 曲線遍歷網(wǎng)格，利用該順序?qū)W(wǎng)格塊進(jìn)行編號和并行分區(qū)，分區(qū)方法利用空間填充Z 曲線的聚類特性，使分在同一個(gè)進(jìn)程的網(wǎng)格塊盡可能保持相鄰，降低了并行發(fā)送與接收方面的需求。