“定黏假設”對伴隨系統求解和梯度精度影響

吳航空，王丁喜，黃秀全，徐慎忍

西北工業大學 動力與能源學院，西安 710072

隨著計算機技術的發展，基于計算流體力學(CFD)的葉輪機設計優化方法獲得了較快發展。目前，得到應用的優化設計方法大致可分為兩類：全局優化方法和局部優化方法。全局優化方法雖然可以獲得全局最優解，但是需要進行百次甚至千次CFD計算。在計算能力及計算資源有限的條件下，其應用受限。相較于全局優化方法，局部優化方法雖然只能獲得局部最優解，但是其在計算耗時方面具有較大優勢。此外，由具有一定經驗研究人員設計的葉輪機葉片往往很接近最優結果，這種局部最優解甚至就是全局最優結果。

一般而言，局部優化方法需要獲得目標函數關于設計變量的梯度信息，也稱目標函數的靈敏度。有限差分方法和線化方法是最早得到應用的兩種靈敏度計算方法。前者的靈敏度計算精度與擾動步長相關。擾動步長太大，截斷誤差會對計算結果精度產生較大影響；擾動步長太小，消去誤差將是誤差的主要來源。后者通過求解線化方程獲得靈敏度信息，其靈敏度計算精度高且不受擾動步長的影響。然而，這兩種方法的共同缺點在于其CFD計算次數與設計變量的維度線性相關。高負荷和高效率等的設計要求使得葉輪機葉片具有復雜的三維彎、扭、掠結構，參數化該葉片外形往往需要成百上千個設計變量，意味著成百上千次CFD計算。對于葉輪機內部多設計變量優化問題，有限差分及線化方法的計算耗時仍然不可接受。相比于上述兩種方法，伴隨方法是一種更加高效的梯度計算方法。每步尋優過程只需分別求解一次非線性流場和線性伴隨場，而與設計變量的個數無關。對于多設計變量優化問題，這大大減少了優化所需要的計算耗時，適合日常工程應用。

伴隨方法有兩種形式：連續伴隨方法和離散伴隨方法,具體的操作流程如圖1所示。連續伴隨方法首先由微分形式的控制方程，經過數學推導獲得微分形式的伴隨方程，然后離散求解。而離散伴隨方法是在離散控制方程的基礎上獲得離散形式的伴隨方程。由于離散伴隨方法的離散格式與離散的流動控制方程具有更好的一致性，靈敏度計算精度高，因此本研究主要基于離散伴隨方法展開。

圖1 連續伴隨方法與離散伴隨方法操作流程圖Fig.1 Operation flow chart of continuous and discrete adjoint methods

離散伴隨程序的開發可以通過手動微分實現，也可以借助自動微分軟件。手動微分獲得的離散伴隨程序效率高，在計算時間和內存消耗方面具有一定優勢，但是手動微分高精度格式和湍流模型方程等較為復雜和困難。此外，手動微分無法自動實現流場求解器和伴隨求解器的同步調整。當對流場求解器的某一子程序進行修改時，需重新推導公式并微分這些子程序，過程繁瑣，容易出錯。自動微分可以完全實現伴隨代碼開發的自動化及流場代碼與伴隨代碼的同步調整，程序的可維護性強。由于整個過程不需要手動微分代碼，程序開發的工作量少，難度低。

隨著自動微分技術的不斷成熟，借助自動微分開發離散伴隨求解器逐漸成為一種趨勢。國外利用自動微分軟件開發離散伴隨求解器的應用較早。Giles等在2000年左右就開始應用自動微分軟件TAPENADE微分流場求解器的子程序，然后再手動組裝這些微分子程序。相比于對整個流場代碼進行微分處理，此方法得到的離散伴隨求解器效率高，計算耗時短。Albring等在開源軟件SU2的基礎上，借助自動微分軟件開發出SU2的伴隨功能，并以復步長變量法的結果作為參考，校驗了該伴隨求解器在梯度計算中的有效性和可靠性。國內，張朝磊等基于自動微分技術構建了離散伴隨優化系統，并對二維無黏跨聲速透平葉柵進行了氣動優化設計。在保證質量流量不變的情況下，優化后出口熵增率減少8.82%，由此驗證了該無黏離散伴隨系統的正確性及有效性。

目前，在以中文形式公開發表的文獻中還未發現基于自動微分技術的全三維湍流伴隨的研究工作。無論是羅佳奇等的連續伴隨方法還是馬燦等的手動離散伴隨方法，都未考慮黏性系數變化對靈敏度精度的影響，即采用“定黏假設”方法。“定黏假設”的引入可簡化伴隨方程的推導及子程序的微分過程，但是會引起靈敏度計算誤差。在某些情況下，這些誤差甚至會改變靈敏度的正負號，使得優化朝著完全錯誤的方向進行。另一方面，目前所研究的“定黏假設”方法同時忽略了層流及湍流黏性系數的影響，因而無法確定層流和湍流黏性系數對靈敏度精度影響程度。在實際應用過程中是否可以選擇性凍結層流或者湍流黏性系數，在保證靈敏度精度的前提下，節省計算時間和存儲。

針對上述問題，本文在自動微分軟件TAPENADE逆向模式的框架下，介紹全三維湍流伴隨求解器的開發過程包括子程序的微分及手動組裝過程，特別是與黏性流動相關子程序的微分與組裝。伴隨靈敏度的驗證采用線化求解器，而線化求解器的開發則是利用自動微分的正向模式。當伴隨求解器與線化求解器都完全微分的時候，兩種方法的雅可比矩陣互為轉置，特征值相同，因此應該具有一致的靈敏度收斂性和相同的漸近收斂率。在此基礎上，本文以跨聲速NASA Rotor 67為研究對象，在兩個不同工況點(近失速點及最高效率點)對比完全湍流伴隨求解器與線化求解器的靈敏度精度、靈敏度收斂性和殘差的漸近收斂率，并且分析研究不同“定黏假設”方法對伴隨求解器計算結果的影響。

1 控制方程

三維圓柱坐標系下雷諾平均Navier-Stokes(RANS)方程為

(1)

式中:

式中:為密度;為壓強;為總內能；為總焓;、、分別為軸向、徑向、周向速度;、、、、、為切應力;、、為正應力;為溫度;為偽時間;、、分別為軸向、徑向、周向熱傳導量；和分別為總的黏性系數和熱傳導系數,它們的值是其層流與湍流分量之和：

(2)

其中：下標l表示層流分量；t表示湍流分量。層流黏性系數由Sutherland定律確定：

(3)

湍流黏性系數通過求解湍流模型方程獲得，本文采用一方程Spalart-Allmaras(SA)湍流模型。SA模型的控制方程為

(4)

(5)

得到黏性系數后，可通過式(6)求解熱傳導系數：

(6)

式中：為定壓比熱比；為普朗特數，其中為0.70，為0.90。

式(1)的時空離散格式為：對流項的離散采用中心差分加標量人工黏性；物理耗散項的離散采用高斯公式將體積分轉化為面積分；偽時間項的離散采用顯隱混合格式，即顯式多階龍格-庫塔和隱式LU-SGS(Lower-Upper Symmetric Gauss-Seidel)。

為了簡化后續推導過程，式(1)可寫成半離散格式：

(7)

式中：為對流項殘差；為黏性項殘差。

2 伴隨方法

葉輪機內流優化問題的數學表達式為

min=(,,)

(8)

s.t.(,,)+(,,,,,)=

(9)

式(8)和式(9)分別為目標函數和約束。為設計變量;為計算域內部流場變量;為邊界處的流場變量;為流場變量的空間一階導數。由于計算域內部和邊界處的流場變量在程序中的處理不同，因此將其分解為和。計算域內部和邊界處的流場變量之間滿足如下顯式關系：

=(,)

(10)

黏性系數與設計變量及流場變量滿足如下關系：

=(,,)=l,t

(11)

流場變量的空間一階導數與設計變量及流場變量之間滿足如下關系：

=(,,)

(12)

線化式(8)可得目標函數關于設計變量的靈敏度：

(13)

式中:和關于的靈敏度可通過求解方程式(9)的線化形式獲得：

(14)

其中:

線化式(10)可得：

(15)

線化式(11)可得：

(16)

線化式(12)可得：

(17)

聯立方程式(14)～式(17)可得：

(18)

式中:

由方程(18)可得：

(19)

將式(15)和式(19)分別代入靈敏度計算式(13) 可得：

(,inv+,vis)(,inv+,vis)

(20)

(21)

線化方程(18)和伴隨方程(21)的共同點在于都為線性方程，并且由于兩種方法的雅可比矩陣互為轉置，特征值相同，因此在線化求解器與伴隨求解器都完全微分的前提下，兩者應具有一致的靈敏度收斂曲線及相同的殘差漸近收斂率。不同的是，線化方程(18)的右邊項與設計變量相關，設計變量的個數決定了所需求解的線化方程的個數。伴隨方程(21)的右邊項與目標函數相關而與設計變量無關，所需求解的伴隨方程個數與目標函數個數線性相關。對于葉輪機內流優化問題，設計變量的維度遠遠大于目標函數的維度，因而伴隨方法更加高效。

在伴隨方法發展的早期，為了簡化推導過程，研究人員通常采用“定黏假設”方法，即不考慮黏性系數(包括層流和湍流)變化對伴隨方程的影響，相應的雅可比矩陣的黏性部分,vis可簡化為

(22)

,vis簡化為

(23)

接下來將分別給出凍結層流黏性系數及湍流黏性系數時的雅可比矩陣。

1) 凍結湍流黏性系數，,vis和,vis為

(24)

(25)

2) 凍結層流黏性系數，,vis和,vis為

(26)

(27)

得到伴隨方程的雅可比矩陣后，可通過時間推進方法迭代求解方程(21)獲得伴隨變量，然后將伴隨變量代入式(20)得到目標函數靈敏度：

(28)

3 流場求解器及伴隨求解器

本節將分別介紹流場求解器和伴隨求解器的數據結構及子程序功能。

3.1 流場求解器

圖2展示了流場求解器的流程圖。要實現其功能，需要6個步驟，每步的具體功能如下：

INIT: 初始化流場。

圖2 非線性流場求解器流程圖Fig.2 Flow chart of nonlinear flow solver

Rinv: 計算對流項殘差。

Viscous residual: 計算黏性項殘差。黏性殘差的計算包含4部分，如圖3所示。

1) Ux: 計算流場變量的空間一階導數。

2) SL: 應用Sutherland定律計算層流黏性系數。

3) SA: 求解湍流模型方程得到湍流黏性系數。

4) Rvis: 計算黏性項殘差。

UPDATE: 求解控制方程(7)，更新內部流場變量。

BC: 施加邊界條件，更新邊界處的流場變量。

圖3 黏性殘差計算Fig.3 Calculation of viscous residual

OBJ: 計算目標函數。

上述6個步驟循環往復直到結果滿足收斂要求。

3.2 湍流伴隨求解器

本研究采用自動微分軟件TAPENADE的逆向模式開發湍流伴隨求解器，具體開發過程為：首先利用自動微分技術對流場求解器的主要子程序進行微分，然后對相應微分子程序按照一定的順序手動組裝。由于逆向模式的原因，手動組裝伴隨求解器的過程較為復雜。圖4展示了湍流伴隨求解器的流程圖，下標a表示由自動微分逆向模式得到的變量或者子程序。從圖中可以看出，湍流伴隨求解器包含兩部分：伴隨方程的求解及靈敏度計算。伴隨方程由于與設計變量的微分無關，因此在這一部分只微分流場變量，凍結設計變量。而靈敏度計算部分，需要微分設計變量。

圖4 湍流伴隨求解器流程圖Fig.4 Flow chart of adjoint solver

3.2.1 伴隨方程求解

伴隨方程的求解需要7個步驟，每步實現的具體功能如下：

OBJ_A1: 計算目標函數關于和的伴隨方向導數，并將其存儲在和，可得：

(29)

式中:為輸入變量，為目標函數的加權因子，只有一個目標函數時，通常將其設為1。

BC_A1: 將轉化到中并累加(方框中+=為累加符號)到，可得：

(30)

對比發現，方程(30)的-與伴隨方程(21)的右邊項一致。

INIT_A: 初始化伴隨變量。為了保證線化求解器與湍流伴隨求解器具有一致的靈敏度收斂曲線，需要將伴隨變量初始化為0。

Rinv_A: 計算伴隨方程對流項殘差。

(31)

Adjoint Viscous Residual: 計算伴隨方程的黏性項殘差。黏性殘差的計算包含4部分，如圖5所示。

圖5 伴隨黏性殘差計算Fig.5 Calculation of adjoint viscous residual

對比圖4和圖5發現，流場求解器的黏性項殘差計算順序與伴隨求解器相反，這是自動微分逆向模式所導致的，因此在組裝與黏性項殘差計算相關微分子程序的時候，需要特別注意。

1) Rvis_A: 計算黏性項殘差，并將其分別存儲在、、、和。

(32)

(33)

2) Ux_A: 將存儲在a的黏性項殘差轉化到和，并累加。

(34)

(35)

3) SL_A: 將存儲在的黏性項殘差轉化到和，并累加。

(36)

(37)

4) SA_A: 將存儲在的黏性項殘差轉化到和，并累加。

(38)

(39)

相加黏性項殘差和對流項殘差得到總殘差：

(40)

BC_A: 將存儲在邊界處的殘差轉化到流場內部并累加。

(41)

對比發現，方程(41)中的與伴隨方程(21) 的左邊項一致。

UPDATE_A: 更新伴隨變量。

3.2.2 靈敏度計算

得到更新的伴隨變量后，將其代入式(27)就可得到目標函數的靈敏度。具體數據結構如圖6 所示。

圖6 靈敏度計算Fig.6 Sensitivity evaluation

OBJ_A2: 計算目標函數關于和的伴隨方向導數，可得：

(42)

注意：該微分子程序只需要被計算一次。將其放在靈敏度計算部分是為了方便讀者理解，在實際程序開發過程中，該微分子程序通常被放在INIT_A(見圖4)微分子程序之前。

Rinv_A2: 計算對流項殘差關于和的伴隨方向導數，可得：

(43)

Adjoint Viscous Residual 2: 計算黏性殘差關于和的伴隨方向導數，此過程包含4步，如圖7所示。

圖7 完全湍流伴隨方法黏性殘差計算Fig.7 Calculation of adjoint viscous residual for ADJ

1) Rvis_A2: 計算黏性殘差關于、、、和的伴隨方向導數。

2) Ux_A2: 計算關于和的伴隨方向導數，并累加。

3) SL_A2: 計算關于和的伴隨方向導數，并累加。

4) SA_A2: 計算關于和的伴隨方向導數，并累加。

上述4個步驟得到的和為

(44)

(45)

相加對應項可得：

(46)

計算靈敏度。

(47)

式中:的轉置為所需求解的目標函數的靈敏度。

3.3 定黏假設

本節將介紹3種不同的“定黏假設”方法：① 凍 結層流及湍流黏性系數(FLEV)；② 凍結層流黏性系數(FLV)；③ 凍結湍流黏性系數(FEV)。為了節省空間，本節只對“定黏假設”伴隨求解器與完全湍流伴隨求解器存在差異的地方展開敘述。

3.3.1 凍結層流和湍流黏性系數

同時凍結層流和湍流黏性系數時，黏性殘差的計算(見圖5)將簡化為圖8。同理，靈敏度部分關于黏性殘差的計算也做相應調整。為了節省空間，不再贅述。

圖8 凍結層流及湍流黏性系數時黏性殘差計算Fig.8 Calculation of adjoint viscous residual with FLEV

3.3.2 凍結層流黏性系數

凍結層流黏性系數時，與黏性殘差計算相關的微分子程序將簡化為圖9。靈敏度計算部分做同步調整。

圖9 凍結層流黏性系數時黏性殘差計算Fig.9 Calculation of adjoint viscous residual with FLV

3.3.3 凍結湍流黏性系數

凍結湍流黏性系數時，與黏性殘差計算相關的子程序將簡化為圖10,同步調整靈敏度計算部分。

圖10 凍結湍流黏性系數時黏性殘差計算Fig.10 Calculation of adjoint viscous residual with FEV

4 流場求解器驗證

本文以跨聲速NASA Rotor 67為研究對象。NASA Rotor 67是NASA Lewis研究中心設計的雙級壓氣機的第一級轉子，其設計參數如表1 所示。

表1 跨聲速NASA Rotor 67的設計參數Table 1 Design parameters of transonic NASA Rotor 67

首先做網格無關性驗證，采用3套不同網格，分別為Grid1，Grid2和Grid3。Grid2的網格如圖11所示，其周向和徑向的網格點數都為65，軸向的網格點數為145(葉片表面的網格點數為81)，總的網格點數大約為613 000。相比于Grid2，Grid1和Grid3的網格點總數分別減小和增加一倍。

進出口為亞聲速邊界條件,在進口給定總溫、總壓及兩個方向的氣流角。總溫及總壓的大小分別為101 325 Pa和288.15 K，氣流角為0°。在出口通過徑向平衡關系獲得背壓沿徑向的分布。為了避免解析黏性剪切底層，在壁面處采用滑移邊界條件及壁面函數。對于葉片前緣之前和尾緣之后的周期性幾何上采用周期性邊界條件。

圖11 NASA Rotor 67算例的計算網格Fig.11 Computational grid of NASA Rotor 67

通過改變出口背壓，進行一系列計算，得到NASA Rotor 67在100%轉速下的特性線。圖12(a) 為流量-壓比特性圖，圖12(b)為流量-效率特性圖。兩幅圖的橫坐標通過實驗測得的堵塞流量(34.96 kg/s)無量綱化得到。從圖12可以看出，Grid1的結果與Grid2和Grid3的結果存在差異，但Grid2和Grid3的計算結果基本重合。為了節省計算時間和計算存儲，接下來的靈敏度驗證部分采用Grid2。

圖12 100%轉速下NASA Rotor 67算例的特性線Fig.12 Performance characteristics of NASA Rotor 67 at 100% design speed

5 結 果

本文選取數值結果的最高效率點及近失速點兩個不同工況(如圖12所示)驗證全三維湍流伴隨求解器的穩定性、靈敏度精度、靈敏度收斂性及殘差的漸近收斂率，并與“定黏假設”結果和線化求解器結果進行對比。最高效率點與近失速點的性能參數如表2所示。

伴隨求解器的驗證部分只考慮性能參數關于邊界條件(入口總溫)的靈敏度。這樣做的好處在于靈敏度計算部分不需要對與幾何變量相關的參數進行微分，簡單、易于操作，并且由于驗證過程不涉及網格擾動，因而可以避免由于網格擾動所造成的誤差。性能參數包括出口質量流量和面積平均總壓。其中，前者為RANS方程守恒變量的線性函數，而后者為RANS方程守恒變量的非線性函數。

表2 NASA Rotor 67最高效率點與近失速點性能參數

5.1 最高效率點

圖13給出了最高效率點98%葉高處的相對馬赫數云圖。從圖中可以看到，通道內有較強激波，并且激波前的馬赫數為1.4，激波的存在使得流場具有較強的非線性。這種復雜非線性流場給湍流伴隨求解器的驗證工作帶來挑戰，但同時也能更進一步驗證湍流伴隨求解器的有效性及魯棒性。

圖13 最高效率點98%葉高的相對馬赫數云圖Fig.13 Relative Mach number contours at 98% span at a peak efficiency point

5.1.1 流場及伴隨場

圖14(a)為50%葉高處的湍流變量云圖，圖14(b) 為與湍流方程所對應的伴隨變量云圖。對比兩幅圖發現，流場變量的下游為伴隨變量的上游，反之亦然。這一特征在一定程度上能用來對湍流伴隨求解器進行定性的校核。

圖14 50%葉高處湍流變量及與湍流方程對應的伴隨變量Fig.14 Flow turbulence variables and adjoint variables corresponding to turbulence model equation at 50% span

5.1.2 靈敏度收斂性及精度

圖15分別展示了質量流量d/d與總壓靈敏度d/d收斂曲線。圖中，ADJ表示完全湍流伴隨求解器，ADJ(FLV)表示凍結層流黏性系數的伴隨求解器，ADJ(FEV)表示凍結湍流黏性系數的伴隨求解器，ADJ(FLEV)表示同時凍結層流和湍流黏性系數的伴隨求解器，LIN表示線化求解器。

以線化求解器結果作為參考，ADJ的靈敏度相對誤差不超過千分之一(如表3所示)。而“定黏假設”方法對伴隨靈敏度精度及收斂性有較大影響。相比于FLV，FEV對伴隨靈敏度的精度影響更大。前者影響不超過1%，而后者的影響超過了5%。對于質量流量靈敏度，FLV會使得伴隨靈敏度增大(相比于LIN)，而FEV會使得伴

圖15 最高效率點伴隨求解器靈敏度與線化求解器靈敏度收斂曲線Fig.15 Sensitivity convergence histories of adjoint and linear solvers near a peak efficiency point

表3 最高效率點伴隨靈敏度相對誤差

隨靈敏度減小，并且影響程度遠遠大于FLV。當采用FLEV的時候，FLV和FEV的影響會相互抵消一部分，因而ADJ(FLEV)的相對誤差反而小于ADJ(FEV)，但是仍然大于ADJ(FLV)。此外，相比于質量流量靈敏度，FLEV對總壓靈敏度的影響更大。采用FLEV，質量流量靈敏度的相對誤差為4.30%，而總壓靈敏度的相對誤差為7.19%。綜上可得，FLEV對非線性目標函數的影響要大于線性目標函數。

5.1.3 殘差的漸近收斂率

在展示殘差收斂曲線之前，首先介紹一下漸近收斂率的定義。式(7)的迭代格式可寫為

+1=-Δτ()

(48)

式中:為總的殘差，包括對流項殘差與黏性項殘差;為迭代步數。假設為控制方程(48)的精確解，則數值結果與精確解的誤差為：

(49)

將方程(49)代入方程(48)，并做一階泰勒展開可得：

(50)

(51)

即矩陣-Δτ的最大特征值決定著誤差的漸近收斂率。

線化方程(18)的誤差漸近收斂率的表達式與非線性方程(51)的完全一致。為了節省空間不再贅述。接下來只針對伴隨方程(21)推導其誤差的漸近收斂率表達式。離散伴隨方程(21)可得：

(52)

假設為控制方程(52)的精確解，則數值結果與精確解之間的誤差為

(53)

因為與互為轉置，特征值相同，因此矩陣-Δ與矩陣-Δτ的特征值相同，最大特征值也相同。由漸近收斂率定義可知，方程(53) 的漸近收斂率與方程(51)相同。從數值結果的角度分析就是完全湍流伴隨殘差與線化殘差的收斂曲線的斜率相同。當采用“定黏假設”方法的時候，由于其影響伴隨方程的雅可比矩陣的特征值，從而影響漸近收斂率，使得伴隨殘差與線化殘差的收斂曲線不一致。

圖16展示了線化殘差與伴隨殘差的質量流量和總壓靈敏殘差收斂曲線(RMS為所有變量的均方根誤差，并取對數)。圖16(b)中的黑色虛線為ADJ殘差收斂曲線的平均斜率曲線。從圖中可以看出，ADJ與LIN的殘差收斂曲線的斜率相同，即一致的漸近收斂率。FLV對漸近收斂率影響較小，ADJ(FLV)的殘差收斂曲線與ADJ重合。而FEV和FLEV對漸近收斂率影響較大。無論是質量流量還是總壓靈敏度，FEV和FLEV都使得殘差收斂曲線斜率的絕對值變小，收斂變慢。

圖16 最高效率點伴隨求解器與線化求解器的殘差收斂曲線Fig.16 Residual convergence histories of adjoint and linear solvers at a peak efficiency point

圖17(a)和圖17(b)分別為50%葉高處的層流黏性系數和湍流黏性系數云圖。從圖中可以看出，層流黏性系數的數值在10量級，并且變化范圍較小(1.6×10～2.0×10)，而湍流黏性系數的數值更大(10量級)，變化范圍也更廣(0.001～0.01)。因而湍流黏性系數對總的黏性系數的貢獻更大，影響也更大，驗證了上述FEV比FLV對伴隨靈敏度影響更大的結論。

圖17 50%葉高層流黏性系數與湍流黏性系數云圖Fig.17 Laminar and eddy viscosity coefficient contours at 50% span

5.1.4 計算時間及計算存儲

圖18對比了不同方法的計算耗時，橫坐標為迭代步數，縱坐標為計算時間。從圖中可以看出，所有伴隨求解器的計算時間和存儲消耗都高于LIN。其中，ADJ(FEV)的計算時間和存儲消耗要低于ADJ(FLV)(見表4)。這是因為湍流黏性系數的計算牽扯到湍流模型方程的求解，計算存儲和計算時間消耗更大，凍結湍流黏性系數將節省更多的計算存儲和計算時間。

圖18 伴隨求解器與線化求解器計算耗時對比Fig.18 Comparison of computational cost between adjoint and linear solvers

表4 伴隨與線化求解器的計算耗時及計算存儲對比

5.2 近失速點

圖19給出了近失速點98%葉高處的相對馬赫數云圖。相較于最高效率點，近失速點的流場更加復雜。不僅有激波的存在，而且大攻角導致吸力面產生較大流動分離。這種同時具有強非線性和流動分離的復雜流場能更進一步驗證湍流伴隨求解器的有效性及魯棒性。

圖19 近失速點98%葉高處的相對馬赫數云圖Fig.19 Relative Mach number contours at 98% span near the stall point

為了節省空間，接下來只針對總壓靈敏度展開分析。近失速點，伴隨求解器與線化求解器的殘差收斂曲線如圖20所示。與設計點相同，ADJ與LIN具有相同的殘差漸近收斂率，并且FLV對計算結果影響較小。不同于設計點，在近失速點采用FEV將直接導致計算結果發散。圖21給出了伴隨求解器與線化求解器的靈敏度收斂曲線。ADJ(FEV)和ADJ(FLEV)的靈敏度都未收斂，而ADJ與LIN的靈敏度收斂曲線吻合度非常高，相對誤差為0.10%，稍高于最高效率點的0.07%。這進一步驗證了湍流伴隨求解器的有效性及魯棒性，也說明了在伴隨求解器開發過程中，考慮湍流黏性系數的微分可以增加伴隨求解器的魯棒性。

圖20 近失速點伴隨求解與線化求解器殘差收斂曲線Fig.20 Residual convergence histories of adjoint and linear solvers at a near stall point

圖21 近失速點伴隨與線化求解器靈敏度收斂曲線Fig.21 Sensitivity convergence histories of adjoint and linear solvers at a near stall point

6 結 論

本文采用自動微分技術開發全三維湍流離散伴隨求解器，并給出相應的流程圖。不僅詳細推導了完全湍流和3種不同定黏方法的伴隨方程，而且介紹了伴隨求解器每一個微分子程序的功能及手動組裝這些微分子程序的過程。最后以NASA Rotor 67為研究對象，分別在最高效率點與近失速點研究了定黏方法對離散伴隨系統的靈敏度精度、靈敏度收斂性及殘差漸近收斂率的影響，并與完全湍流伴隨求解器及線化求解器的結果進行對比，得到的結論如下：

1) ADJ與LIN具有一致的靈敏度收斂曲線及相同殘差漸近收斂率。無論是線性目標函數還是非線性目標函數，相對誤差都不超過0.10%。

2) “定黏假設”方法會影響伴隨求解器的靈敏度精度及收斂性，并且FEV的影響大于FLV，前者對靈敏度的影響超過5%，而后者的影響低于1%。

3) “定黏假設”方法會影響伴隨求解器的殘差漸近收斂率，并且FEV的影響遠大于FLV。在最高效率點，FEV和FLEV會使得伴隨求解器的殘差收斂變慢；在近失速點，FEV和FLEV將直接導致伴隨求解器的計算結果發散。

4) 相比于ADJ，由“定黏假設”方法得到的伴隨求解器在計算時間和計算存儲消耗方面具有一定的優勢，但是仍然遠高于LIN。

致 謝

感謝課題組張倩為本文提供NASA Rotor 67的幾何數據及實驗數據。