旋轉(zhuǎn)坐標(biāo)系下分區(qū)計(jì)算的LU隱式方法

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

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

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

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

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

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

1 數(shù)值計(jì)算方法

1.1 LU-SGS時(shí)間隱式迭代格式

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

(1)

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

λij(Qj-Qi))sij

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

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

(2)

式中：

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

得到

同時(shí)用到

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

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

1.2 基于LU分解的時(shí)間迭代并行計(jì)算方法

由1.1節(jié)可見，高斯-賽德爾迭代方法在并行分區(qū)邊界會(huì)碰到問題，無(wú)法做到與串行計(jì)算結(jié)果一致。而在本質(zhì)上，DP-LUR方法將離散方程式(2) 的非對(duì)角項(xiàng)移到方程的右邊，用Jacobian迭代代替對(duì)稱Gauss-Seidel迭代。這樣DP-LUR方法對(duì)當(dāng)前時(shí)間層的流場(chǎng)值沒有依賴，保證了計(jì)算域內(nèi)部單元和邊界單元采用完全相同的計(jì)算格式推進(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

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

ΔQ0=D-1R

Dok=1,2,…,ksweep

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

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

End Do

ΔQn+1=ΔQ(ksweep)

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

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

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

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

1.3 并行程序設(shè)計(jì)

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

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

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

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

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

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

……

上述數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)中只考慮了計(jì)算單元，沒有考慮虛擬邊界單元。當(dāng)計(jì)算每一個(gè)具體的網(wǎng)格塊時(shí)，另開辟專門的邊界數(shù)組進(jìn)行存儲(chǔ)虛擬單元的流場(chǎng)變量值。如圖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 計(jì)算結(jié)果與分析

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

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

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

并行計(jì)算采用8個(gè)分區(qū)，網(wǎng)格劃分方式分為兩種，一種是沿展向劃分，一種是沿周向劃分(見圖6)。為了考察分區(qū)邊界對(duì)數(shù)值計(jì)算的影響，分別對(duì)3種時(shí)間離散方法進(jìn)行了對(duì)比，一種是LU-SGS 方法，每個(gè)時(shí)間步分別進(jìn)行一次L掃描和U掃描；一種是DP-LUR方法，該方法由于迭代效率較低，在每個(gè)時(shí)間步的內(nèi)迭代步數(shù)為8步；一種是HLU-SGS方法，該方法結(jié)合了LU-SGS和DP-LUR方法，每個(gè)時(shí)間步進(jìn)行2次DP-LUR迭代和1次LU-SGS迭代。需要進(jìn)一步指出的是，由于DP-LUR方法采用雅克比內(nèi)迭代，文獻(xiàn)[8]認(rèn)為5次內(nèi)迭代可以獲得與高斯-賽德爾迭代類似的收斂速度，本文測(cè)試結(jié)果表明，DP-LUR方法采用8次內(nèi)迭代，HLU-SGS方法采用2次DP-LUR迭代和1次LU-SGS迭代可以獲得很高穩(wěn)定性，使周向邊界取20倍旋翼半徑、CFL=10 000情況下保持計(jì)算穩(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種隱式時(shí)間迭代方法得到的殘差變化曲線，其中殘差用密度的L2模表示，Span8表示展向分為8個(gè)分區(qū)，Circum8表示周向分為8個(gè)分區(qū)。從圖中可以看到，對(duì)于DP-LUR格式而言，無(wú)論網(wǎng)格是周向剖分還是展向剖分，密度殘差一直保持為機(jī)器零值，即網(wǎng)格運(yùn)動(dòng)沒有引起流場(chǎng)發(fā)生變化，流場(chǎng)保持為靜止初場(chǎng)。對(duì)于LU-SGS格式來(lái)說(shuō)，展向剖分網(wǎng)格對(duì)應(yīng)的流場(chǎng)隨迭代步數(shù)不變，但是對(duì)于周向剖分網(wǎng)格，隨迭代步數(shù)密度殘差增大，流場(chǎng)不再是理論上的靜止流場(chǎng)。為什么是周向劃分網(wǎng)格導(dǎo)致計(jì)算發(fā)散呢？對(duì)于當(dāng)前旋轉(zhuǎn)運(yùn)動(dòng)網(wǎng)格來(lái)說(shuō)，周向剖分網(wǎng)格的邊界方向與旋轉(zhuǎn)方向相同，邊界面上的逆變速度(流體相對(duì)于網(wǎng)格坐標(biāo)的速度)大，也就意味著面上的流通量大。當(dāng)LU-SGS格式用于上述網(wǎng)格的并行計(jì)算時(shí)，由于其在邊界上的簡(jiǎn)化處理，方程中起主導(dǎo)作用的物理量被丟掉，導(dǎo)致了誤差放大和計(jì)算發(fā)散。相對(duì)而言，HLU-SGS和DP-LUR格式在時(shí)間迭代中增加了邊界上的信息交換，對(duì)兩種剖分網(wǎng)格的計(jì)算都得到了收斂結(jié)果。

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

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

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

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

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

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

為了測(cè)試當(dāng)前的計(jì)算方法，選擇標(biāo)準(zhǔn)模型Caradonna-Tung旋翼[17]進(jìn)行了懸停狀態(tài)計(jì)算。該旋翼包括兩片等弦長(zhǎng)無(wú)扭轉(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個(gè)方向的維數(shù)為109×119×165(對(duì)應(yīng)徑向、周向和法向)，槳尖區(qū)域的最小網(wǎng)格間距為0.05倍弦長(zhǎng)，網(wǎng)格單元數(shù)約200萬(wàn)；密網(wǎng)格3個(gè)方向的網(wǎng)格維數(shù)為305×177×969(對(duì)應(yīng)徑向、周向和法向)，如圖10所示，將槳葉展向0.65R～1.1R范圍并且槳葉下方區(qū)域作為槳尖渦模擬區(qū)對(duì)網(wǎng)格進(jìn)行均勻加密，網(wǎng)格間距為0.02倍弦長(zhǎng)，網(wǎng)格單元數(shù)超過5千萬(wàn)。計(jì)算采用256個(gè)分區(qū)的并行計(jì)算，為了達(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

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

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

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

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

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

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

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

圖17給出了稀/密兩套網(wǎng)格的截面流線和下洗速度云圖比較，可以看到兩套網(wǎng)格得到的下洗場(chǎng)分布比較一致，最大下洗速度約為15 m/s。不同之處是密網(wǎng)格刻畫了更多的流動(dòng)細(xì)節(jié)，包括旋翼尾跡中槳尖渦的分布，以及由于槳尖渦的存在使局部流線發(fā)生卷曲等。另外，值得注意的是，密網(wǎng)格計(jì)算得到的第一個(gè)槳尖渦的強(qiáng)度大于稀網(wǎng)格，旋渦引起的局部上洗速度大于15 m/s，由于旋渦位置靠近槳葉，將對(duì)槳葉氣動(dòng)性能的預(yù)測(cè)產(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給出的是槳尖渦運(yùn)動(dòng)軌跡的計(jì)算值與試驗(yàn)值比較。從圖中可以看到，密網(wǎng)格模擬的槳尖渦下降軌跡與試驗(yàn)值吻合一致，稀網(wǎng)格則過高估計(jì)了槳尖渦下降速度；從槳尖渦向內(nèi)收縮來(lái)看，計(jì)算與Caradonna和Tung的試驗(yàn)值[17]存在較大差異。原因分析如下，本文計(jì)算只模擬了旋翼槳葉，不包括槳轂和旋翼軸；試驗(yàn)數(shù)據(jù)取自1981年Caradonna和Tung的試驗(yàn)，試驗(yàn)裝置不僅有旋翼軸，還有測(cè)壓管線沿旋翼軸引出，使旋翼軸體偏大，目測(cè)半徑有5%～10%倍旋翼半徑，中心軸體部分會(huì)影響槳尖渦向內(nèi)收縮的偏移量，定性地說(shuō)中心軸體增大，槳尖渦向內(nèi)收縮量減小。圖18同時(shí)參考了Kocurek的尾跡模型數(shù)據(jù)[17]，與本文計(jì)算結(jié)果比較一致，一定程度上驗(yàn)證了本文計(jì)算?？偟膩?lái)說(shuō)，當(dāng)前密網(wǎng)格計(jì)算獲得了旋翼懸停流場(chǎng)的主要流動(dòng)特征，在旋翼氣動(dòng)性能預(yù)測(cè)上更有優(yōu)勢(shì)，通過采用適用于分區(qū)計(jì)算的LU隱式方法改進(jìn)以后，稀網(wǎng)格和密網(wǎng)格都獲得了充分收斂流場(chǎng)，進(jìn)一步驗(yàn)證了本文的并行分區(qū)算法。

圖18 槳尖渦的運(yùn)動(dòng)軌跡Fig.18 Trajectories of blade tip vortex motion

3 結(jié) 論

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

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

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

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