肖中云，劉剛，牟斌，江雄

中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

當(dāng)前CFD并行計算大多采用區(qū)域分解的并行算法，將計算域按照負(fù)載平衡原則分配給不同進(jìn)程，每個進(jìn)程分別采用時間推進(jìn)方法進(jìn)行迭代求解，并通過進(jìn)程間消息傳遞交換流場數(shù)據(jù)。這種方法由于信息交換量小、負(fù)載平衡性好而得到廣泛使用。嚴(yán)格地說，隱式時間迭代算法由于在分區(qū)邊界存在依賴關(guān)系并不能直接用于并行，通常這種影響量很小而被忽略。但是這樣的簡化在旋翼計算中卻碰到了收斂性變差的問題，甚至導(dǎo)致一些非物理的計算結(jié)果。出現(xiàn)這樣的問題與旋翼流動特點有關(guān)，旋翼以旋轉(zhuǎn)運動為主，在懸停和低速前飛狀態(tài)下，旋翼尾跡旋繞槳盤附近，向遠(yuǎn)場的發(fā)展過程緩慢。由于流動的這種特性，數(shù)值計算誤差容易被積累，導(dǎo)致收斂性變差[1-3]。隨著并行計算規(guī)模的增長，旋翼CFD計算對算法的可擴(kuò)展性提出了很高要求，這就需要發(fā)展考慮分區(qū)邊界影響的隱式時間迭代格式，對忽略邊界影響的處理方法進(jìn)行重新評估，目的是使大規(guī)模并行計算不會因為進(jìn)程數(shù)增加導(dǎo)致收斂性變差或者失效，同時計算結(jié)果不受進(jìn)程數(shù)目的影響。

文獻(xiàn)[4-5]對計算流體力學(xué)的并行計算技術(shù)進(jìn)行了綜述，指出顯式格式是內(nèi)在并行的而且并行網(wǎng)格劃分不影響其數(shù)值特性，保持串、并行計算過程一致。常見隱式格式中雅克比松弛是內(nèi)在并行的，高斯-賽德爾松弛與其他線松弛格式由于網(wǎng)格點相互依賴，非內(nèi)在并行，但通常有更好的收斂性。由于這類方法的并行求解比較困難，實際計算中通常采取“局部隱”的方法，即對計算區(qū)域進(jìn)行分塊，對每塊獨立進(jìn)行隱式求解，塊邊界值由相鄰塊前一時刻的值來確定。

基于LU分解的對稱高斯-賽德爾(LU-SGS)迭代是一種典型的隱式時間迭代算法[6-7]。Wissink等[8]最早提出了原始LU-SGS時間迭代方法在分區(qū)邊界會產(chǎn)生誤差，尤其在并行規(guī)模較分區(qū)邊界較多時誤差得到放大；提出了基于雅克比迭代的數(shù)據(jù)并行LU松弛算法(DP-LUR)，同時將LU-SGS與DP-LUR進(jìn)行結(jié)合構(gòu)造了混合LU-SGS方法(HLU-SGS)，提高了收斂效率。Luo等[9]在對HLU-SGS方法進(jìn)行研究以后，認(rèn)為當(dāng)網(wǎng)格塊數(shù)較多時，混合方法具有優(yōu)勢，當(dāng)并行規(guī)模較小時邊界數(shù)據(jù)傳遞不會對結(jié)果帶來影響，由于LU-SGS具有較高的迭代效率，在計算中直接采用了LU-SGS方法進(jìn)行并行計算。Otero和Eliasson[10]將LU-SGS直接應(yīng)用于并行計算，采用多次迭代的方法，在每一次前向掃描和后向掃描完成后分別進(jìn)行一次消息傳遞接口(MPI)交換信息，更新邊界單元的ΔQ值。文中稱這種方法為混沌更新，利用多次迭代消除未按順序更新帶來的誤差。

目前在旋翼計算中隱式時間迭代方法得到了廣泛應(yīng)用，但是關(guān)于分區(qū)邊界影響的研究較少。文獻(xiàn)[11]在非慣性旋轉(zhuǎn)坐標(biāo)系下對旋翼懸停狀態(tài)進(jìn)行了模擬，時間離散采用的LU-SGS隱式迭代方法，在兩次掃描過程中，虛擬單元上的值直接賦零。文獻(xiàn)[12-13]在慣性系下對旋翼前飛狀態(tài)進(jìn)行并行計算，時間推進(jìn)采用隱式LU-SGS格式，但對并行計算的邊界如何處理沒有進(jìn)行特別說明。針對旋翼并行計算中隱式時間迭代方面的潛在問題，本文在同一框架下發(fā)展了LU-SGS、DP-LUR和HLU-SGS計算方法，分析了原始LU-SGS方法在并行計算中產(chǎn)生誤差的行為及原因，通過在旋轉(zhuǎn)坐標(biāo)系下靜止場計算驗證了算法的有效性。在此基礎(chǔ)上，采用分區(qū)并行計算、多重網(wǎng)格方法等先進(jìn)算法對旋翼的懸停流場進(jìn)行模擬，獲得了較高的并行效率和流場計算精度，進(jìn)一步驗證了本文方法對于大規(guī)模分區(qū)計算的適用性。相關(guān)方法可以進(jìn)一步推廣到其他流動問題的旋轉(zhuǎn)坐標(biāo)系計算。

1 數(shù)值計算方法

1.1 LU-SGS時間隱式迭代格式

Navier-Stokes方程經(jīng)過隱式時間離散后的形式為[6]

(1)

式中：Q為守恒型流場變量；R為Navier-Stokes方程右端項；V為單元控制體體積；Δt為時間步長；下標(biāo)i為單元編號；上標(biāo)n為當(dāng)前時間步。為了得到左端的雅克比矩陣，通量函數(shù)采用最大特征值分裂，即

λij(Qj-Qi))sij

式中：F為通矢量；j為單元i的相鄰單元；sij為交界面面積；nij為單位法向矢量；λij為通量雅可比矩陣的譜半徑，λij=|vij·nij|+cij。

將左端雅克比矩陣分解為上三角矩陣U、下三角矩陣L和對角陣D這3個部分，得到

(2)

式中：

其中：N(i)為單元i的相鄰單元；L(i)為左相鄰單元；U(i)為右相鄰側(cè)單元；由于控制體封閉[14]，單元i滿足：

得到

同時用到

LU-SGS算法首先對式(2)進(jìn)行近似分解，然后采用高斯-賽德爾迭代求解線性代數(shù)方程，即

(D+L)D-1(D+U)ΔQ=R

1) L掃描

(D+L)ΔQ*=R

2) U掃描

(D+U)ΔQ=DΔQ*

式中：ΔQ*為L掃描后得到的流場變量值；ΔQn+1為U掃描后得到的流場變量值。MPI并行計算采用域分解的算法，高斯-賽德爾迭代在域邊界位置會碰到相鄰點流場值無法及時更新的問題。圖1 給出的是二維網(wǎng)格LU-SGS迭代的L掃描示意圖。當(dāng)掃描到邊界單元時，位于其左側(cè)或下側(cè)相鄰單元的當(dāng)前時刻值(ΔQ*)不可用，通常只能進(jìn)行簡化處理，如假設(shè)ΔQ*=0，這種簡化帶來的誤差在特定情況下可能放大，引起計算不穩(wěn)定或者發(fā)散。在并行計算分區(qū)數(shù)很大時，對邊界的處理方式要謹(jǐn)慎對待。

圖1 LU-SGS格式L掃描示意圖Fig.1 Diagram of lower sweep in LU-SGS scheme

1.2 基于LU分解的時間迭代并行計算方法

由1.1節(jié)可見，高斯-賽德爾迭代方法在并行分區(qū)邊界會碰到問題，無法做到與串行計算結(jié)果一致。而在本質(zhì)上，DP-LUR方法將離散方程式(2) 的非對角項移到方程的右邊，用Jacobian迭代代替對稱Gauss-Seidel迭代。這樣DP-LUR方法對當(dāng)前時間層的流場值沒有依賴，保證了計算域內(nèi)部單元和邊界單元采用完全相同的計算格式推進(jìn)。DP-LUR算法的收斂效率低于LU-SGS算法，通常將迭代次數(shù)取為5步甚至更多。

ΔQ0=D-1R

Dok=1,2,…,ksweep

ΔQk=D-1[R-(U+L)ΔQk-1]

End Do

ΔQn+1=ΔQksweep

為了同時利用LU-SGS迭代的高效率和DP-LUR在邊界處理上的優(yōu)點，文獻(xiàn)[8]發(fā)展了兩者相混合的方法(HLU-SGS)，即在塊邊界單元采用Jacobian迭代格式，計算邊界單元的ΔQ，隨后對所有采用LU-SGS方法進(jìn)行計算，對于塊內(nèi)部單元來說，算法等同于LU-SGS方法；對于塊邊界單元來說，相鄰單元的ΔQ為已知的當(dāng)前時刻的值。在極限情況下，當(dāng)所有單元都是邊界單元的情況時，該計算格式退化為點隱式的Jacobian迭代方法。HLU-SGS的算法為

ΔQ0=D-1R

Dok=1,2,…,ksweep

邊界信息交換ΔQ*(k)=ΔQ(k-1)

采用Gauss-Seidel方法計算ΔQ(k)

End Do

ΔQn+1=ΔQ(ksweep)

圖2給出的是多塊結(jié)構(gòu)網(wǎng)格的計算流程圖，其中MGCycle表示多重網(wǎng)格循環(huán)，BLK表示網(wǎng)格塊號，P1、P2、PM表示進(jìn)程編號?？梢钥吹矫窟M(jìn)行一次時間步迭代，各進(jìn)程間只進(jìn)行一次通信，交換網(wǎng)格邊界的的流場信息。在每個計算進(jìn)程內(nèi)，各網(wǎng)格塊之間相對獨立，為了節(jié)省內(nèi)存開銷，大部分工作變量都定義在當(dāng)前網(wǎng)格塊上，包括流場變量、網(wǎng)格度量系數(shù)、右端項殘差等，在完成當(dāng)前塊計算以后，這些數(shù)據(jù)被下一網(wǎng)格塊取代。

并行隱式時間迭代算法的流程圖如圖3所示，從圖中可以看到，各網(wǎng)格塊之間的計算相互耦合，比如在每次進(jìn)行時間迭代時，需要對所有網(wǎng)格塊進(jìn)行掃描，并進(jìn)行進(jìn)程間通信，交換網(wǎng)格邊界流場解ΔQ信息。顯然，在上述需求下，原來的基于單塊的數(shù)組存取方式不再適合。本文的改進(jìn)方式如下：取消基于單塊網(wǎng)格的數(shù)組變量，用一維數(shù)組存放當(dāng)前進(jìn)程下所有網(wǎng)格的流場變量；開辟全局?jǐn)?shù)組存放網(wǎng)格的幾何度量系數(shù)信息，這些信息在右端項計算和時間迭代計算中將要用到；開辟全局?jǐn)?shù)組存放右端項殘差，該信息在DP-LUR迭代和LU-SGS 迭代中將分別用到。在經(jīng)過上述改進(jìn)以后，在同一計算框架下，既可以完成原有的LU-SGS迭代計算，也可以進(jìn)行多次內(nèi)迭代的DP-LUR計算，還可以進(jìn)行DP-LUR和LU-SGS的混合計算。

圖2 網(wǎng)格分區(qū)并行計算的原有流程Fig.2 Original flow chart of grid partitioned parallel computation

圖3 網(wǎng)格分區(qū)并行計算的改進(jìn)流程Fig.3 Improved flow chart of grid partitioned parallel computation

1.3 并行程序設(shè)計

并行網(wǎng)格分配采用貪婪優(yōu)化算法，其核心思想是：將當(dāng)前最重的任務(wù)，分配給當(dāng)前負(fù)載最小的進(jìn)程[15]。當(dāng)原始網(wǎng)格分布不滿足分配要求時，需要對網(wǎng)格塊進(jìn)行剖分，根據(jù)貪婪算法的設(shè)計思想，在找到當(dāng)前最需要網(wǎng)格的進(jìn)程和當(dāng)前最大網(wǎng)格塊之后，如果網(wǎng)格塊單元數(shù)小于或剛好等于要求的網(wǎng)格單元數(shù)時，則直接將該網(wǎng)格塊放到該進(jìn)程。但更多時候是，當(dāng)前網(wǎng)格塊單元數(shù)大于所需的網(wǎng)格單元數(shù)，就需要對當(dāng)前網(wǎng)格塊進(jìn)行拆分。采取的拆分策略是：當(dāng)網(wǎng)格塊的剖分位置位于維數(shù)25%～75%之間時，對網(wǎng)格進(jìn)行按需剖分并分配；當(dāng)剖分位置不在25%～75%范圍內(nèi)時，對網(wǎng)格塊進(jìn)行中位剖分，剖分后不分配，繼續(xù)查找當(dāng)前最大網(wǎng)格塊進(jìn)行下一步操作。負(fù)載平衡與邊界單元數(shù)最少是結(jié)構(gòu)網(wǎng)格剖分的一對矛盾，為了使網(wǎng)格塊不被分得過于零碎，同時也減少分區(qū)邊界數(shù)，要犧牲一些負(fù)載平衡要求，一般將負(fù)載不平衡度控制在5%～7%以內(nèi)可以取得較好的剖分結(jié)果。并行網(wǎng)格分配的另一項重要任務(wù)是建立分區(qū)邊界的映射關(guān)系，確定出發(fā)送/接收次數(shù)，以及每次發(fā)送/接收的對方進(jìn)程號、字段長度、相對應(yīng)的網(wǎng)格單元編號等信息。

圖4 MPI并行的編程模型Fig.4 Model for MPI parallel programming

并行程序設(shè)計采用對等模式，如圖4所示，每個進(jìn)程地位相同，代碼和功能基本一致，只不過是處理的數(shù)據(jù)對象不同。由于邊界點新值的計算需要相鄰邊界其他塊的數(shù)據(jù)，因此在每一個數(shù)據(jù)塊的兩側(cè)各增加1列數(shù)據(jù)空間，用于存放從相鄰數(shù)據(jù)塊通信得到的數(shù)據(jù)，進(jìn)程之間通過MPI發(fā)送與接收函數(shù)進(jìn)行消息傳遞。

為了適應(yīng)多重網(wǎng)格計算，網(wǎng)格單元按照原始網(wǎng)格、第2重網(wǎng)格、第3重網(wǎng)格等的順序進(jìn)行排序，并存儲在一維數(shù)組當(dāng)中，每一重網(wǎng)格計算分別是對數(shù)組中的某一段進(jìn)行操作。假設(shè)當(dāng)前進(jìn)程的原始網(wǎng)格單元為N，那么多重網(wǎng)格的次密層單元數(shù)為N/8，第3重為N/64，以此類推。對任意網(wǎng)格點來說，由多重網(wǎng)格層數(shù)、網(wǎng)格塊號、3個方向的指標(biāo)號(i,j,k)可以確定出唯一指標(biāo)號，即在一維數(shù)組中的存放位置。假設(shè)任意流場變量A，在進(jìn)行多重網(wǎng)格計算時，每一重網(wǎng)格上的數(shù)據(jù)分別為

第1重網(wǎng)格：A(1:N)

第2重網(wǎng)格：A(N+1:N+N/8)

……

上述數(shù)據(jù)存儲結(jié)構(gòu)中只考慮了計算單元，沒有考慮虛擬邊界單元。當(dāng)計算每一個具體的網(wǎng)格塊時，另開辟專門的邊界數(shù)組進(jìn)行存儲虛擬單元的流場變量值。如圖5所示，“i”方向的邊界數(shù)組可以表示為“Boundi(1:2,jdim,kdim)”，其中“1，2”分別代表左邊界的兩排虛擬單元，“3,4”代表右邊界的兩排虛擬單元，jdim和kdim表示邊界面j和k方向上的網(wǎng)格維數(shù)。

圖5 網(wǎng)格塊與虛擬邊界單元Fig.5 Grid block and virtual boundary cell

2 計算結(jié)果與分析

2.1 旋轉(zhuǎn)網(wǎng)格算例測試

根據(jù)動網(wǎng)格計算對“均勻流保持特性”的要求[16]，設(shè)計了如下算例，即取靜止流場中的一塊網(wǎng)格，采用動網(wǎng)格方法對流場進(jìn)行模擬，無論網(wǎng)格做何種運動，最終得到的結(jié)果都應(yīng)該是保持靜止初場不變。如圖6所示，取計算域為周向角等于36°的一塊扇形網(wǎng)格，假設(shè)參考長度(旋翼展長)為R，網(wǎng)格的周向遠(yuǎn)場邊界為4R。采用靜止流場作為初場，經(jīng)過任意長時間的計算后，流場應(yīng)保持初場不變，即流場狀態(tài)變量與環(huán)境變量相等，用密度或其他變量表示的流場殘值應(yīng)該保持為機(jī)器零值。

計算在旋轉(zhuǎn)坐標(biāo)系下進(jìn)行，旋翼槳尖馬赫數(shù)Matip=0.352。由于計算域為沒有旋翼的空流場，遠(yuǎn)場為靜止大氣，邊界條件采用無反射黎曼邊界條件，采用靜止流場作為初場，空間離散采用Roe通量差分格式，庫朗數(shù)(CFL數(shù))取為20。首先對單塊網(wǎng)格的串行計算進(jìn)行了測試，結(jié)果表明，在當(dāng)前參數(shù)設(shè)置下計算可以順利進(jìn)行，密度殘差一直保持在機(jī)器零值(實數(shù)用雙精度表示，機(jī)器零值約為10-16)。

并行計算采用8個分區(qū)，網(wǎng)格劃分方式分為兩種，一種是沿展向劃分，一種是沿周向劃分(見圖6)。為了考察分區(qū)邊界對數(shù)值計算的影響，分別對3種時間離散方法進(jìn)行了對比，一種是LU-SGS 方法，每個時間步分別進(jìn)行一次L掃描和U掃描；一種是DP-LUR方法，該方法由于迭代效率較低，在每個時間步的內(nèi)迭代步數(shù)為8步；一種是HLU-SGS方法，該方法結(jié)合了LU-SGS和DP-LUR方法，每個時間步進(jìn)行2次DP-LUR迭代和1次LU-SGS迭代。需要進(jìn)一步指出的是，由于DP-LUR方法采用雅克比內(nèi)迭代，文獻(xiàn)[8]認(rèn)為5次內(nèi)迭代可以獲得與高斯-賽德爾迭代類似的收斂速度，本文測試結(jié)果表明，DP-LUR方法采用8次內(nèi)迭代，HLU-SGS方法采用2次DP-LUR迭代和1次LU-SGS迭代可以獲得很高穩(wěn)定性，使周向邊界取20倍旋翼半徑、CFL=10 000情況下保持計算穩(wěn)定。綜合上述因素，建議DP-LUR方法內(nèi)迭代次數(shù)選取為5～8次，HLU-SGS方法的DP-LUR內(nèi)迭代次數(shù)為1～2次。

圖6 展向和周向網(wǎng)格剖分示意圖Fig.6 Schematics of spanwise and circumferential grid partitions

圖7給出的是2種剖分網(wǎng)格、3種隱式時間迭代方法得到的殘差變化曲線，其中殘差用密度的L2模表示，Span8表示展向分為8個分區(qū)，Circum8表示周向分為8個分區(qū)。從圖中可以看到，對于DP-LUR格式而言，無論網(wǎng)格是周向剖分還是展向剖分，密度殘差一直保持為機(jī)器零值，即網(wǎng)格運動沒有引起流場發(fā)生變化，流場保持為靜止初場。對于LU-SGS格式來說，展向剖分網(wǎng)格對應(yīng)的流場隨迭代步數(shù)不變，但是對于周向剖分網(wǎng)格，隨迭代步數(shù)密度殘差增大，流場不再是理論上的靜止流場。為什么是周向劃分網(wǎng)格導(dǎo)致計算發(fā)散呢？對于當(dāng)前旋轉(zhuǎn)運動網(wǎng)格來說，周向剖分網(wǎng)格的邊界方向與旋轉(zhuǎn)方向相同，邊界面上的逆變速度(流體相對于網(wǎng)格坐標(biāo)的速度)大，也就意味著面上的流通量大。當(dāng)LU-SGS格式用于上述網(wǎng)格的并行計算時，由于其在邊界上的簡化處理，方程中起主導(dǎo)作用的物理量被丟掉，導(dǎo)致了誤差放大和計算發(fā)散。相對而言，HLU-SGS和DP-LUR格式在時間迭代中增加了邊界上的信息交換，對兩種剖分網(wǎng)格的計算都得到了收斂結(jié)果。

圖7 3種格式、2種網(wǎng)格下的殘差歷程(CFL=20)Fig.7 Residual histories of three operators applied to two grids (CFL=20)

針對圖7中的發(fā)散算例(對應(yīng)周向剖分網(wǎng)格)，圖8為采用不同CFL數(shù)計算得到的殘差曲線。從圖中可以看到，在迭代2 000步前提下，LU-SGS格式CFL=20對應(yīng)的密度殘差發(fā)散最快，CFL=10對應(yīng)的曲線次之，CFL=5對應(yīng)的曲線在計算時間內(nèi)仍一直保持收斂狀態(tài)，后續(xù)計算會不會發(fā)散未知。由此可見，這種由計算格式帶入的誤差十分危險，可能是以一種不易察覺的方式使計算結(jié)果受到污染。對于HLU-SGS格式來說，每個時間步采用2次DP-LUR迭代和1次LU-SGS迭代的組合方式使計算格式非常魯棒，在CFL=10 000情況下計算仍保持收斂。

從前面分析可以看到，造成計算不穩(wěn)定的原因來自于網(wǎng)格分區(qū)邊界的運動速度，在旋轉(zhuǎn)坐標(biāo)系下該最大速度取決于轉(zhuǎn)速和遠(yuǎn)場邊界的大小。下面測試保持旋轉(zhuǎn)速度不變，將遠(yuǎn)場邊界與旋轉(zhuǎn)中心距離b從4R分別增大到10R和20R，對應(yīng)周向網(wǎng)格運動速度為1.41、3.52和7.04倍聲速。圖9給出的是不同計算域大小對收斂性的影響，其中時間推進(jìn)為HLU-SGS格式，CFL=10 000。從圖中可以看到，在遠(yuǎn)場邊界不斷增大的情況下，HLU-SGS計算格式都達(dá)到了收斂狀態(tài)，顯示了良好的計算穩(wěn)定性；另外隨著遠(yuǎn)場邊界的增加，殘差量級增大，在遠(yuǎn)場邊界為20R時殘差只能收斂到10-14量級，說明旋轉(zhuǎn)坐標(biāo)系計算遠(yuǎn)場邊界不宜取過大。

圖8 CFL數(shù)對殘差曲線的影響Fig.8 Effects of CFL number on residual histories

圖9 計算域?qū)埐钋€的影響(HLU-SGS)Fig.9 Effects of computational domain on residual histories (HLU-SGS)

2.2 Caradonna-Tung標(biāo)模旋翼計算

為了測試當(dāng)前的計算方法，選擇標(biāo)準(zhǔn)模型Caradonna-Tung旋翼[17]進(jìn)行了懸停狀態(tài)計算。該旋翼包括兩片等弦長無扭轉(zhuǎn)的槳葉，翼剖面形狀為NACA-0012翼型，槳盤半徑為1.143 m，槳葉展弦比為6。模擬狀態(tài)槳尖馬赫數(shù)Matip=0.44，總距θc=8°。網(wǎng)格采用重疊網(wǎng)格方法，包括槳葉網(wǎng)格和背景網(wǎng)格，背景網(wǎng)格又分別生成了稀/密程度不同的兩套網(wǎng)格，其中稀網(wǎng)格3個方向的維數(shù)為109×119×165(對應(yīng)徑向、周向和法向)，槳尖區(qū)域的最小網(wǎng)格間距為0.05倍弦長，網(wǎng)格單元數(shù)約200萬；密網(wǎng)格3個方向的網(wǎng)格維數(shù)為305×177×969(對應(yīng)徑向、周向和法向)，如圖10所示，將槳葉展向0.65R～1.1R范圍并且槳葉下方區(qū)域作為槳尖渦模擬區(qū)對網(wǎng)格進(jìn)行均勻加密，網(wǎng)格間距為0.02倍弦長，網(wǎng)格單元數(shù)超過5千萬。計算采用256個分區(qū)的并行計算，為了達(dá)到負(fù)載平衡要求，網(wǎng)格經(jīng)過了多次剖分，剖分后的稀網(wǎng)格為452塊，密網(wǎng)格為448塊，圖10顯示了背景網(wǎng)格在經(jīng)過網(wǎng)格剖分后產(chǎn)生的分區(qū)邊界。

圖10 稀/密網(wǎng)格的并行分區(qū)Fig.10 Grid partitions of coarse and fine grids

計算參考系為原點位于旋翼中心的旋轉(zhuǎn)坐標(biāo)系，只需模擬單片槳葉，周向采用周期性邊界條件，遠(yuǎn)場邊界采用懸停源/匯邊界條件[18-19]。流動假設(shè)為全湍流，湍流模型為兩方程k-ω剪切應(yīng)力輸運(SST)模型，空間離散采用Roe通量差分格式。時間迭代分別采用了原始LU-SGS方法、數(shù)據(jù)并行的DP-LUR方法與混合方法HLU-SGS，其中DP-LUR方法每個時間步進(jìn)行6次內(nèi)迭代，HLU-SGS方法每個時間步進(jìn)行1次DP-LUR迭代和1次LU-SGS迭代。前面提到原始LU-SGS方法有計算不穩(wěn)定的情況發(fā)生，為了對比3種計算方法，這里選擇采用較小計算域和低CFL數(shù)方法，使3種隱式格式都能夠正常計算，具體為背景網(wǎng)格的周向邊界距離旋翼軸3倍旋翼半徑，CFL=20。此外，為了加速收斂，計算還采用了V循環(huán)的多重網(wǎng)格方法和低速預(yù)處理方法，計算目標(biāo)是使流場達(dá)到充分收斂。

流場收斂性研究與格式比較在稀網(wǎng)格上進(jìn)行，為了充分讓流場收斂，總迭代步數(shù)達(dá)到10萬步。圖11給出的是拉力系數(shù)CT/σ隨迭代步數(shù)的變化，其中CT表示用槳盤面積無量綱化的拉力系數(shù)，σ表示旋翼實度，可以看到3種隱式算法中DP-LUR收斂最快，其他兩種方法相當(dāng)，HLU-SGS略快于LU-SGS方法。由于DP-LUR方法在一個時間步內(nèi)進(jìn)行了6次迭代，其收斂快是可以預(yù)料的，HLU-SGS方法在一個時間步數(shù)內(nèi)進(jìn)行了1次DP-LUR迭代和1次LU-SGS迭代，收斂略快于原有的LU-SGS方法。

由于3種方法每個時間步的計算時間不同，為了比較各自的收斂速度，圖12給出了拉力系數(shù)CT/σ隨時間的變化曲線，可以看到原始的LU-SGS方法收斂速度最快，DP-LUR方法收斂速度最慢，盡管HLU-SGS方法相對原始LU-SGS方法計算量有增加，這是邊界隱式計算所付出的代價，但是相對DP-LUR方法計算效率有了明顯提高。

圖11 拉力系數(shù)隨迭代步數(shù)的變化Fig.11 Variations of thrust coefficient with iterations

圖12 拉力系數(shù)隨時間的變化Fig.12 Variations of thrust coefficient with time

圖13給出的是并行計算的加速比和并行效率統(tǒng)計，在高性能計算機(jī)系統(tǒng)上完成計算，每個計算節(jié)點采用Intel Xeon處理器(64核，256 G共享內(nèi)存)，結(jié)點之間采用無限帶寬技術(shù)(Infiniband)高速互連。測試算例為密網(wǎng)格算例(網(wǎng)格單元數(shù)約5千萬)，在同樣計算條件下分別采用了1～512核進(jìn)行了計算，由于所測得的數(shù)據(jù)有一定的散布，圖中除了給出每個算例的測試結(jié)果外，還用曲線進(jìn)行了擬合。從圖中可以看到，在測試范圍內(nèi)，隨并行規(guī)模的擴(kuò)大，加速比單調(diào)增大，但并行效率呈下降趨勢。以128核并行計算為例，加速比約為94.5倍，并行效率約為74%。由于本文采用了并行隱式時間迭代方法，相對原始LU-SGS方法，在每個時間步增加了一次信息交換(傳遞ΔQ值)，通信需求的增加使并行效率有一定下降。

圖13 并行計算的加速比與效率Fig.13 Speedup ratio and efficiency of parallel computations

圖17給出了稀/密兩套網(wǎng)格的截面流線和下洗速度云圖比較，可以看到兩套網(wǎng)格得到的下洗場分布比較一致，最大下洗速度約為15 m/s。不同之處是密網(wǎng)格刻畫了更多的流動細(xì)節(jié)，包括旋翼尾跡中槳尖渦的分布，以及由于槳尖渦的存在使局部流線發(fā)生卷曲等。另外，值得注意的是，密網(wǎng)格計算得到的第一個槳尖渦的強度大于稀網(wǎng)格，旋渦引起的局部上洗速度大于15 m/s，由于旋渦位置靠近槳葉，將對槳葉氣動性能的預(yù)測產(chǎn)生影響。

圖14 拉力與扭矩系數(shù)的收斂曲線Fig.14 Convergence histories of thrust and torque coefficients

圖15 Q值等值面Fig.15 Iso-surface of Q value

圖16 稀/密網(wǎng)格的渦量云圖比較Fig.16 Comparison of vorticity contour of coarse and fine grids

圖17 截面流線和下洗速度云圖Fig.17 Surface streamlines and down wash velocity contour

圖18給出的是槳尖渦運動軌跡的計算值與試驗值比較。從圖中可以看到，密網(wǎng)格模擬的槳尖渦下降軌跡與試驗值吻合一致，稀網(wǎng)格則過高估計了槳尖渦下降速度；從槳尖渦向內(nèi)收縮來看，計算與Caradonna和Tung的試驗值[17]存在較大差異。原因分析如下，本文計算只模擬了旋翼槳葉，不包括槳轂和旋翼軸；試驗數(shù)據(jù)取自1981年Caradonna和Tung的試驗，試驗裝置不僅有旋翼軸，還有測壓管線沿旋翼軸引出，使旋翼軸體偏大，目測半徑有5%～10%倍旋翼半徑，中心軸體部分會影響槳尖渦向內(nèi)收縮的偏移量，定性地說中心軸體增大，槳尖渦向內(nèi)收縮量減小。圖18同時參考了Kocurek的尾跡模型數(shù)據(jù)[17]，與本文計算結(jié)果比較一致，一定程度上驗證了本文計算?？偟膩碚f，當(dāng)前密網(wǎng)格計算獲得了旋翼懸停流場的主要流動特征，在旋翼氣動性能預(yù)測上更有優(yōu)勢，通過采用適用于分區(qū)計算的LU隱式方法改進(jìn)以后，稀網(wǎng)格和密網(wǎng)格都獲得了充分收斂流場，進(jìn)一步驗證了本文的并行分區(qū)算法。

圖18 槳尖渦的運動軌跡Fig.18 Trajectories of blade tip vortex motion

3 結(jié) 論

本文采用3種基于LU分解的隱式時間迭代算法，基于網(wǎng)格分區(qū)和基于MPI的并行方法，在旋轉(zhuǎn)坐標(biāo)系下對單獨旋轉(zhuǎn)網(wǎng)格、旋翼的懸停狀態(tài)進(jìn)行了數(shù)值模擬，研究了不同時間迭代格式對旋翼計算收斂性的影響。

1) LU-SGS隱式時間迭代方法由于在分區(qū)邊界缺少相鄰單元的流場信息，理論上如果直接應(yīng)用于并行計算會在邊界上帶來誤差，一般認(rèn)為這種誤差量很小、對結(jié)果沒有顯著影響。本文通過旋轉(zhuǎn)網(wǎng)格算例證明，將LU-SGS方法應(yīng)用于旋轉(zhuǎn)網(wǎng)格并行計算時，當(dāng)分區(qū)邊界與旋轉(zhuǎn)方向同向時，時間離散產(chǎn)生的誤差會起主導(dǎo)作用并不斷積累放大，使計算結(jié)果受到污染，并導(dǎo)致計算發(fā)散。

2) DP-LUR方法用雅克比迭代代替高斯-賽德爾迭代算法，避免了分區(qū)邊界的迭代誤差問題，能消除并行剖分對計算的影響，不足之處是每步時間推進(jìn)需要進(jìn)行多次迭代，計算效率低于LU-SGS方法。HLU-SGS方法是DP-LUR與LU-SGS兩種方法的結(jié)合，能有效消除分區(qū)邊界的影響，同時獲得較高的計算效率，是一種適用于并行的隱式時間推進(jìn)方法。算例表明該方法能很好適用于旋轉(zhuǎn)坐標(biāo)系下的旋翼計算，能夠獲得氣動力與尾跡流場都充分收斂的計算結(jié)果。

3) 本文討論的隱式時間迭代方法適用于所用旋轉(zhuǎn)坐標(biāo)系計算，解決了結(jié)構(gòu)網(wǎng)格多塊分區(qū)和并行計算分區(qū)的問題，適合于計算域較大和網(wǎng)格整體旋轉(zhuǎn)的情況。除旋翼懸停計算以外，該方法還可以推廣到飛機(jī)的繞速度矢滾轉(zhuǎn)、螺旋槳、風(fēng)力機(jī)等其他旋轉(zhuǎn)類流場的數(shù)值計算。