“定黏假設(shè)”對(duì)伴隨系統(tǒng)求解和梯度精度影響

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

西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院，西安 710072

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

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

伴隨方法有兩種形式：連續(xù)伴隨方法和離散伴隨方法,具體的操作流程如圖1所示。連續(xù)伴隨方法首先由微分形式的控制方程，經(jīng)過(guò)數(shù)學(xué)推導(dǎo)獲得微分形式的伴隨方程，然后離散求解。而離散伴隨方法是在離散控制方程的基礎(chǔ)上獲得離散形式的伴隨方程。由于離散伴隨方法的離散格式與離散的流動(dòng)控制方程具有更好的一致性，靈敏度計(jì)算精度高，因此本研究主要基于離散伴隨方法展開(kāi)。

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

離散伴隨程序的開(kāi)發(fā)可以通過(guò)手動(dòng)微分實(shí)現(xiàn)，也可以借助自動(dòng)微分軟件。手動(dòng)微分獲得的離散伴隨程序效率高，在計(jì)算時(shí)間和內(nèi)存消耗方面具有一定優(yōu)勢(shì)，但是手動(dòng)微分高精度格式和湍流模型方程等較為復(fù)雜和困難。此外，手動(dòng)微分無(wú)法自動(dòng)實(shí)現(xiàn)流場(chǎng)求解器和伴隨求解器的同步調(diào)整。當(dāng)對(duì)流場(chǎng)求解器的某一子程序進(jìn)行修改時(shí)，需重新推導(dǎo)公式并微分這些子程序，過(guò)程繁瑣，容易出錯(cuò)。自動(dòng)微分可以完全實(shí)現(xiàn)伴隨代碼開(kāi)發(fā)的自動(dòng)化及流場(chǎng)代碼與伴隨代碼的同步調(diào)整，程序的可維護(hù)性強(qiáng)。由于整個(gè)過(guò)程不需要手動(dòng)微分代碼，程序開(kāi)發(fā)的工作量少，難度低。

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

目前，在以中文形式公開(kāi)發(fā)表的文獻(xiàn)中還未發(fā)現(xiàn)基于自動(dòng)微分技術(shù)的全三維湍流伴隨的研究工作。無(wú)論是羅佳奇等的連續(xù)伴隨方法還是馬燦等的手動(dòng)離散伴隨方法，都未考慮黏性系數(shù)變化對(duì)靈敏度精度的影響，即采用“定黏假設(shè)”方法。“定黏假設(shè)”的引入可簡(jiǎn)化伴隨方程的推導(dǎo)及子程序的微分過(guò)程，但是會(huì)引起靈敏度計(jì)算誤差。在某些情況下，這些誤差甚至?xí)淖冹`敏度的正負(fù)號(hào)，使得優(yōu)化朝著完全錯(cuò)誤的方向進(jìn)行。另一方面，目前所研究的“定黏假設(shè)”方法同時(shí)忽略了層流及湍流黏性系數(shù)的影響，因而無(wú)法確定層流和湍流黏性系數(shù)對(duì)靈敏度精度影響程度。在實(shí)際應(yīng)用過(guò)程中是否可以選擇性凍結(jié)層流或者湍流黏性系數(shù)，在保證靈敏度精度的前提下，節(jié)省計(jì)算時(shí)間和存儲(chǔ)。

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

1 控制方程

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

(1)

式中:

式中:為密度;為壓強(qiáng);為總內(nèi)能；為總焓;、、分別為軸向、徑向、周向速度;、、、、、為切應(yīng)力;、、為正應(yīng)力;為溫度;為偽時(shí)間;、、分別為軸向、徑向、周向熱傳導(dǎo)量；和分別為總的黏性系數(shù)和熱傳導(dǎo)系數(shù),它們的值是其層流與湍流分量之和：

(2)

其中：下標(biāo)l表示層流分量；t表示湍流分量。層流黏性系數(shù)由Sutherland定律確定：

(3)

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

(4)

(5)

得到黏性系數(shù)后，可通過(guò)式(6)求解熱傳導(dǎo)系數(shù)：

(6)

式中：為定壓比熱比；為普朗特?cái)?shù)，其中為0.70，為0.90。

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

為了簡(jiǎn)化后續(xù)推導(dǎo)過(guò)程，式(1)可寫(xiě)成半離散格式：

(7)

式中：為對(duì)流項(xiàng)殘差；為黏性項(xiàng)殘差。

2 伴隨方法

葉輪機(jī)內(nèi)流優(yōu)化問(wèn)題的數(shù)學(xué)表達(dá)式為

min=(,,)

(8)

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

(9)

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

=(,)

(10)

黏性系數(shù)與設(shè)計(jì)變量及流場(chǎng)變量滿足如下關(guān)系：

=(,,)=l,t

(11)

流場(chǎng)變量的空間一階導(dǎo)數(shù)與設(shè)計(jì)變量及流場(chǎng)變量之間滿足如下關(guān)系：

=(,,)

(12)

線化式(8)可得目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的靈敏度：

(13)

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

(14)

其中:

線化式(10)可得：

(15)

線化式(11)可得：

(16)

線化式(12)可得：

(17)

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

(18)

式中:

由方程(18)可得：

(19)

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

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

(20)

(21)

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

在伴隨方法發(fā)展的早期，為了簡(jiǎn)化推導(dǎo)過(guò)程，研究人員通常采用“定黏假設(shè)”方法，即不考慮黏性系數(shù)(包括層流和湍流)變化對(duì)伴隨方程的影響，相應(yīng)的雅可比矩陣的黏性部分,vis可簡(jiǎn)化為

(22)

,vis簡(jiǎn)化為

(23)

接下來(lái)將分別給出凍結(jié)層流黏性系數(shù)及湍流黏性系數(shù)時(shí)的雅可比矩陣。

1) 凍結(jié)湍流黏性系數(shù)，,vis和,vis為

(24)

(25)

2) 凍結(jié)層流黏性系數(shù)，,vis和,vis為

(26)

(27)

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

(28)

3 流場(chǎng)求解器及伴隨求解器

本節(jié)將分別介紹流場(chǎng)求解器和伴隨求解器的數(shù)據(jù)結(jié)構(gòu)及子程序功能。

3.1 流場(chǎng)求解器

圖2展示了流場(chǎng)求解器的流程圖。要實(shí)現(xiàn)其功能，需要6個(gè)步驟，每步的具體功能如下：

INIT: 初始化流場(chǎng)。

圖2 非線性流場(chǎng)求解器流程圖Fig.2 Flow chart of nonlinear flow solver

Rinv: 計(jì)算對(duì)流項(xiàng)殘差。

Viscous residual: 計(jì)算黏性項(xiàng)殘差。黏性殘差的計(jì)算包含4部分，如圖3所示。

1) Ux: 計(jì)算流場(chǎng)變量的空間一階導(dǎo)數(shù)。

2) SL: 應(yīng)用Sutherland定律計(jì)算層流黏性系數(shù)。

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

4) Rvis: 計(jì)算黏性項(xiàng)殘差。

UPDATE: 求解控制方程(7)，更新內(nèi)部流場(chǎng)變量。

BC: 施加邊界條件，更新邊界處的流場(chǎng)變量。

圖3 黏性殘差計(jì)算Fig.3 Calculation of viscous residual

OBJ: 計(jì)算目標(biāo)函數(shù)。

上述6個(gè)步驟循環(huán)往復(fù)直到結(jié)果滿足收斂要求。

3.2 湍流伴隨求解器

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

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

3.2.1 伴隨方程求解

伴隨方程的求解需要7個(gè)步驟，每步實(shí)現(xiàn)的具體功能如下：

OBJ_A1: 計(jì)算目標(biāo)函數(shù)關(guān)于和的伴隨方向?qū)?shù)，并將其存儲(chǔ)在和，可得：

(29)

式中:為輸入變量，為目標(biāo)函數(shù)的加權(quán)因子，只有一個(gè)目標(biāo)函數(shù)時(shí)，通常將其設(shè)為1。

BC_A1: 將轉(zhuǎn)化到中并累加(方框中+=為累加符號(hào))到，可得：

(30)

對(duì)比發(fā)現(xiàn)，方程(30)的-與伴隨方程(21)的右邊項(xiàng)一致。

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

Rinv_A: 計(jì)算伴隨方程對(duì)流項(xiàng)殘差。

(31)

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

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

對(duì)比圖4和圖5發(fā)現(xiàn)，流場(chǎng)求解器的黏性項(xiàng)殘差計(jì)算順序與伴隨求解器相反，這是自動(dòng)微分逆向模式所導(dǎo)致的，因此在組裝與黏性項(xiàng)殘差計(jì)算相關(guān)微分子程序的時(shí)候，需要特別注意。

1) Rvis_A: 計(jì)算黏性項(xiàng)殘差，并將其分別存儲(chǔ)在、、、和。

(32)

(33)

2) Ux_A: 將存儲(chǔ)在a的黏性項(xiàng)殘差轉(zhuǎn)化到和，并累加。

(34)

(35)

3) SL_A: 將存儲(chǔ)在的黏性項(xiàng)殘差轉(zhuǎn)化到和，并累加。

(36)

(37)

4) SA_A: 將存儲(chǔ)在的黏性項(xiàng)殘差轉(zhuǎn)化到和，并累加。

(38)

(39)

相加黏性項(xiàng)殘差和對(duì)流項(xiàng)殘差得到總殘差：

(40)

BC_A: 將存儲(chǔ)在邊界處的殘差轉(zhuǎn)化到流場(chǎng)內(nèi)部并累加。

(41)

對(duì)比發(fā)現(xiàn)，方程(41)中的與伴隨方程(21) 的左邊項(xiàng)一致。

UPDATE_A: 更新伴隨變量。

3.2.2 靈敏度計(jì)算

得到更新的伴隨變量后，將其代入式(27)就可得到目標(biāo)函數(shù)的靈敏度。具體數(shù)據(jù)結(jié)構(gòu)如圖6 所示。

圖6 靈敏度計(jì)算Fig.6 Sensitivity evaluation

OBJ_A2: 計(jì)算目標(biāo)函數(shù)關(guān)于和的伴隨方向?qū)?shù)，可得：

(42)

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

Rinv_A2: 計(jì)算對(duì)流項(xiàng)殘差關(guān)于和的伴隨方向?qū)?shù)，可得：

(43)

Adjoint Viscous Residual 2: 計(jì)算黏性殘差關(guān)于和的伴隨方向?qū)?shù)，此過(guò)程包含4步，如圖7所示。

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

1) Rvis_A2: 計(jì)算黏性殘差關(guān)于、、、和的伴隨方向?qū)?shù)。

2) Ux_A2: 計(jì)算關(guān)于和的伴隨方向?qū)?shù)，并累加。

3) SL_A2: 計(jì)算關(guān)于和的伴隨方向?qū)?shù)，并累加。

4) SA_A2: 計(jì)算關(guān)于和的伴隨方向?qū)?shù)，并累加。

上述4個(gè)步驟得到的和為

(44)

(45)

相加對(duì)應(yīng)項(xiàng)可得：

(46)

計(jì)算靈敏度。

(47)

式中:的轉(zhuǎn)置為所需求解的目標(biāo)函數(shù)的靈敏度。

3.3 定黏假設(shè)

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

3.3.1 凍結(jié)層流和湍流黏性系數(shù)

同時(shí)凍結(jié)層流和湍流黏性系數(shù)時(shí)，黏性殘差的計(jì)算(見(jiàn)圖5)將簡(jiǎn)化為圖8。同理，靈敏度部分關(guān)于黏性殘差的計(jì)算也做相應(yīng)調(diào)整。為了節(jié)省空間，不再贅述。

圖8 凍結(jié)層流及湍流黏性系數(shù)時(shí)黏性殘差計(jì)算Fig.8 Calculation of adjoint viscous residual with FLEV

3.3.2 凍結(jié)層流黏性系數(shù)

凍結(jié)層流黏性系數(shù)時(shí)，與黏性殘差計(jì)算相關(guān)的微分子程序?qū)⒑?jiǎn)化為圖9。靈敏度計(jì)算部分做同步調(diào)整。

圖9 凍結(jié)層流黏性系數(shù)時(shí)黏性殘差計(jì)算Fig.9 Calculation of adjoint viscous residual with FLV

3.3.3 凍結(jié)湍流黏性系數(shù)

凍結(jié)湍流黏性系數(shù)時(shí)，與黏性殘差計(jì)算相關(guān)的子程序?qū)⒑?jiǎn)化為圖10,同步調(diào)整靈敏度計(jì)算部分。

圖10 凍結(jié)湍流黏性系數(shù)時(shí)黏性殘差計(jì)算Fig.10 Calculation of adjoint viscous residual with FEV

4 流場(chǎng)求解器驗(yàn)證

本文以跨聲速NASA Rotor 67為研究對(duì)象。NASA Rotor 67是NASA Lewis研究中心設(shè)計(jì)的雙級(jí)壓氣機(jī)的第一級(jí)轉(zhuǎn)子，其設(shè)計(jì)參數(shù)如表1 所示。

表1 跨聲速NASA Rotor 67的設(shè)計(jì)參數(shù)Table 1 Design parameters of transonic NASA Rotor 67

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

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

圖11 NASA Rotor 67算例的計(jì)算網(wǎng)格Fig.11 Computational grid of NASA Rotor 67

通過(guò)改變出口背壓，進(jìn)行一系列計(jì)算，得到NASA Rotor 67在100%轉(zhuǎn)速下的特性線。圖12(a) 為流量-壓比特性圖，圖12(b)為流量-效率特性圖。兩幅圖的橫坐標(biāo)通過(guò)實(shí)驗(yàn)測(cè)得的堵塞流量(34.96 kg/s)無(wú)量綱化得到。從圖12可以看出，Grid1的結(jié)果與Grid2和Grid3的結(jié)果存在差異，但Grid2和Grid3的計(jì)算結(jié)果基本重合。為了節(jié)省計(jì)算時(shí)間和計(jì)算存儲(chǔ)，接下來(lái)的靈敏度驗(yàn)證部分采用Grid2。

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

5 結(jié) 果

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

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

表2 NASA Rotor 67最高效率點(diǎn)與近失速點(diǎn)性能參數(shù)

5.1 最高效率點(diǎn)

圖13給出了最高效率點(diǎn)98%葉高處的相對(duì)馬赫數(shù)云圖。從圖中可以看到，通道內(nèi)有較強(qiáng)激波，并且激波前的馬赫數(shù)為1.4，激波的存在使得流場(chǎng)具有較強(qiáng)的非線性。這種復(fù)雜非線性流場(chǎng)給湍流伴隨求解器的驗(yàn)證工作帶來(lái)挑戰(zhàn)，但同時(shí)也能更進(jìn)一步驗(yàn)證湍流伴隨求解器的有效性及魯棒性。

圖13 最高效率點(diǎn)98%葉高的相對(duì)馬赫數(shù)云圖Fig.13 Relative Mach number contours at 98% span at a peak efficiency point

5.1.1 流場(chǎng)及伴隨場(chǎng)

圖14(a)為50%葉高處的湍流變量云圖，圖14(b) 為與湍流方程所對(duì)應(yīng)的伴隨變量云圖。對(duì)比兩幅圖發(fā)現(xiàn)，流場(chǎng)變量的下游為伴隨變量的上游，反之亦然。這一特征在一定程度上能用來(lái)對(duì)湍流伴隨求解器進(jìn)行定性的校核。

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

5.1.2 靈敏度收斂性及精度

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

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

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

表3 最高效率點(diǎn)伴隨靈敏度相對(duì)誤差

隨靈敏度減小，并且影響程度遠(yuǎn)遠(yuǎn)大于FLV。當(dāng)采用FLEV的時(shí)候，F(xiàn)LV和FEV的影響會(huì)相互抵消一部分，因而ADJ(FLEV)的相對(duì)誤差反而小于ADJ(FEV)，但是仍然大于ADJ(FLV)。此外，相比于質(zhì)量流量靈敏度，F(xiàn)LEV對(duì)總壓靈敏度的影響更大。采用FLEV，質(zhì)量流量靈敏度的相對(duì)誤差為4.30%，而總壓靈敏度的相對(duì)誤差為7.19%。綜上可得，F(xiàn)LEV對(duì)非線性目標(biāo)函數(shù)的影響要大于線性目標(biāo)函數(shù)。

5.1.3 殘差的漸近收斂率

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

+1=-Δτ()

(48)

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

(49)

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

(50)

(51)

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

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

(52)

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

(53)

因?yàn)?span id="j5i0abt0b" class="emphasis_italic">與互為轉(zhuǎn)置，特征值相同，因此矩陣-Δ與矩陣-Δτ的特征值相同，最大特征值也相同。由漸近收斂率定義可知，方程(53) 的漸近收斂率與方程(51)相同。從數(shù)值結(jié)果的角度分析就是完全湍流伴隨殘差與線化殘差的收斂曲線的斜率相同。當(dāng)采用“定黏假設(shè)”方法的時(shí)候，由于其影響伴隨方程的雅可比矩陣的特征值，從而影響漸近收斂率，使得伴隨殘差與線化殘差的收斂曲線不一致。

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

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

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

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

5.1.4 計(jì)算時(shí)間及計(jì)算存儲(chǔ)

圖18對(duì)比了不同方法的計(jì)算耗時(shí)，橫坐標(biāo)為迭代步數(shù)，縱坐標(biāo)為計(jì)算時(shí)間。從圖中可以看出，所有伴隨求解器的計(jì)算時(shí)間和存儲(chǔ)消耗都高于LIN。其中，ADJ(FEV)的計(jì)算時(shí)間和存儲(chǔ)消耗要低于ADJ(FLV)(見(jiàn)表4)。這是因?yàn)橥牧黟ば韵禂?shù)的計(jì)算牽扯到湍流模型方程的求解，計(jì)算存儲(chǔ)和計(jì)算時(shí)間消耗更大，凍結(jié)湍流黏性系數(shù)將節(jié)省更多的計(jì)算存儲(chǔ)和計(jì)算時(shí)間。

圖18 伴隨求解器與線化求解器計(jì)算耗時(shí)對(duì)比Fig.18 Comparison of computational cost between adjoint and linear solvers

表4 伴隨與線化求解器的計(jì)算耗時(shí)及計(jì)算存儲(chǔ)對(duì)比

5.2 近失速點(diǎn)

圖19給出了近失速點(diǎn)98%葉高處的相對(duì)馬赫數(shù)云圖。相較于最高效率點(diǎn)，近失速點(diǎn)的流場(chǎng)更加復(fù)雜。不僅有激波的存在，而且大攻角導(dǎo)致吸力面產(chǎn)生較大流動(dòng)分離。這種同時(shí)具有強(qiáng)非線性和流動(dòng)分離的復(fù)雜流場(chǎng)能更進(jìn)一步驗(yàn)證湍流伴隨求解器的有效性及魯棒性。

圖19 近失速點(diǎn)98%葉高處的相對(duì)馬赫數(shù)云圖Fig.19 Relative Mach number contours at 98% span near the stall point

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

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

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

6 結(jié) 論

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

1) ADJ與LIN具有一致的靈敏度收斂曲線及相同殘差漸近收斂率。無(wú)論是線性目標(biāo)函數(shù)還是非線性目標(biāo)函數(shù)，相對(duì)誤差都不超過(guò)0.10%。

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

3) “定黏假設(shè)”方法會(huì)影響伴隨求解器的殘差漸近收斂率，并且FEV的影響遠(yuǎn)大于FLV。在最高效率點(diǎn)，F(xiàn)EV和FLEV會(huì)使得伴隨求解器的殘差收斂變慢；在近失速點(diǎn)，F(xiàn)EV和FLEV將直接導(dǎo)致伴隨求解器的計(jì)算結(jié)果發(fā)散。

4) 相比于ADJ，由“定黏假設(shè)”方法得到的伴隨求解器在計(jì)算時(shí)間和計(jì)算存儲(chǔ)消耗方面具有一定的優(yōu)勢(shì)，但是仍然遠(yuǎn)高于LIN。

致 謝

感謝課題組張倩為本文提供NASA Rotor 67的幾何數(shù)據(jù)及實(shí)驗(yàn)數(shù)據(jù)。