黃江濤, 劉 剛, 周 鑄,*, 高正紅, 黃 勇(. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽 6000; . 西北工業(yè)大學(xué), 陜西 西安 7007)

基于離散伴隨方程求解梯度信息的若干問題研究

黃江濤1, 劉 剛1, 周 鑄1,*, 高正紅2, 黃 勇1
(1. 中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽 621000; 2. 西北工業(yè)大學(xué), 陜西 西安 710072)

基于自主研發(fā)的大規(guī)模并行化結(jié)構(gòu)化網(wǎng)格RANS求解器PMB3D，開展了黏性離散伴隨方程構(gòu)造、求解方法的研究與討論。首先對(duì)離散伴隨求解梯度的思想進(jìn)行簡(jiǎn)要介紹，進(jìn)一步對(duì)無黏項(xiàng)、人工黏性項(xiàng)、黏性項(xiàng)部分對(duì)離散伴隨方程貢獻(xiàn)以及變分推導(dǎo)進(jìn)行了詳細(xì)介紹；文中對(duì)離散伴隨方程無黏項(xiàng)、黏性項(xiàng)邊界條件實(shí)現(xiàn)形式進(jìn)行了詳細(xì)研究，并對(duì)關(guān)鍵模塊變分推導(dǎo)的一些簡(jiǎn)化方式進(jìn)行了研究討論，通過典型寬體飛機(jī)標(biāo)模、外壓式超聲速進(jìn)氣道算例，分析了所采用的簡(jiǎn)化處理方式對(duì)不同問題梯度求解精度的影響。 最后在并行化求解、時(shí)間推進(jìn)以及加速收斂方面進(jìn)行了探討、驗(yàn)證。數(shù)值模擬表明，文中采用的離散伴隨方程形式更有利于程序化、模塊化，梯度計(jì)算精度完全滿足氣動(dòng)優(yōu)化設(shè)計(jì)需要。

離散伴隨方程；氣動(dòng)優(yōu)化；梯度信息；邊界條件；時(shí)間推進(jìn)；并行計(jì)算

0 引 言

基于流動(dòng)變分思想的分析手段以其獨(dú)有的優(yōu)勢(shì)，在氣動(dòng)設(shè)計(jì)、網(wǎng)格誤差修正領(lǐng)域等扮演著重要角色，針對(duì)不同形式的主控方程，CFD學(xué)者們發(fā)展出了連續(xù)伴隨、離散伴隨方程[1-5]，并進(jìn)行了相應(yīng)的求解方式研究，以期高效地獲得氣動(dòng)特性對(duì)設(shè)計(jì)變量的梯度信息，由于該項(xiàng)技術(shù)求解梯度信息的工作量幾乎與設(shè)計(jì)變量個(gè)數(shù)無關(guān)，因此，倍受CFD研究人員以及氣動(dòng)優(yōu)化設(shè)計(jì)研究人員的重視。其中，由于離散伴隨方程與NS方程清晰的導(dǎo)數(shù)關(guān)系，實(shí)現(xiàn)起來比較方便，梯度信息更為準(zhǔn)確等優(yōu)點(diǎn)，在伴隨系統(tǒng)研究領(lǐng)域最為關(guān)注，也是國內(nèi)外空氣動(dòng)力學(xué)研究機(jī)構(gòu)重點(diǎn)發(fā)展的研究方向，大多研究機(jī)構(gòu)均基于自身研發(fā)的大型并行CFD計(jì)算代碼發(fā)展了離散伴隨優(yōu)化系統(tǒng)，例如NASA Langley研究中心利用自動(dòng)微分工具開發(fā)了基于結(jié)構(gòu)化求解器CFL3D、非結(jié)構(gòu)化求解器FUN3D的離散伴隨優(yōu)化系統(tǒng)[6-7]；德宇航基于結(jié)構(gòu)化求解器Flower、非結(jié)構(gòu)化求解器TAU發(fā)展了離散伴隨優(yōu)化系統(tǒng)[8]，法宇航基于CFD代碼elsA開發(fā)了離散伴隨優(yōu)化[9],英國謝菲爾德大學(xué)覃寧[10]開展了基于結(jié)構(gòu)化網(wǎng)格的并行離散伴隨優(yōu)化。

國內(nèi)在離散伴隨方程求解器自主研發(fā)方面也取得了一定的進(jìn)展，例如西北工業(yè)大學(xué)左英桃、高正紅基于結(jié)構(gòu)化網(wǎng)格求解器開展了M6機(jī)翼離散伴隨優(yōu)化[11]；熊俊濤[12]等基于顯式時(shí)間推進(jìn)實(shí)現(xiàn)了離散伴隨方程的求解；屈崑等利用Tapenade自動(dòng)微分工具進(jìn)行通量變分，按照矩陣模式組裝到全局稀疏矩陣，實(shí)現(xiàn)了穩(wěn)態(tài)CFD的伴隨系統(tǒng)求解[13]；西安交通大學(xué)張朝磊等基于離散伴隨理論和自動(dòng)微分技術(shù)構(gòu)建離散伴隨系統(tǒng)，應(yīng)用于透平葉柵的氣動(dòng)優(yōu)化[14]；南京航空航天大學(xué)高宜勝、伍貽兆、夏健等基于非結(jié)構(gòu)求解器進(jìn)行了翼型離散伴隨優(yōu)化[15]；中國空氣動(dòng)力研究與發(fā)展中心李彬等基于非結(jié)構(gòu)求解器實(shí)現(xiàn)了離散伴隨系統(tǒng)的開發(fā)[16]。

對(duì)于離散伴隨優(yōu)化設(shè)計(jì)系統(tǒng)的研究，大部分研究側(cè)重點(diǎn)處于氣動(dòng)優(yōu)化問題本身，而對(duì)于離散伴隨系統(tǒng)構(gòu)造的關(guān)鍵問題歸納整理較少，本文將系統(tǒng)地對(duì)離散伴隨方程求解梯度信息的若干問題開展研究；另一方面，離散伴隨優(yōu)化多點(diǎn)設(shè)計(jì)的方式是加權(quán)平均，設(shè)計(jì)結(jié)果依賴于權(quán)系數(shù)的選擇，需要對(duì)權(quán)系數(shù)進(jìn)行多次嘗試，由于結(jié)構(gòu)化網(wǎng)格求解器在相同資源條件下計(jì)算效率、精度更高，更利于變權(quán)系數(shù)的多點(diǎn)分布式優(yōu)化設(shè)計(jì)，比較適合于工程型號(hào)設(shè)計(jì)。因此，本文基于中國空氣動(dòng)力研究與發(fā)展中心自主研發(fā)的大型并行結(jié)構(gòu)化網(wǎng)格RANS求解器PMB3D[17]，開展了離散伴隨系統(tǒng)構(gòu)造的關(guān)鍵問題研究。

1 基于伴隨思想的梯度計(jì)算方法簡(jiǎn)介

對(duì)于氣動(dòng)優(yōu)化設(shè)計(jì)問題：

及其殘差約束R(W,X,D)=0，可以構(gòu)造以下目標(biāo)函數(shù)：

對(duì)式(2)進(jìn)行求導(dǎo)，

式(4)就是流場(chǎng)伴隨方程，通過迭代方法求解Λ之后，可以通過式(6)進(jìn)行梯度信息快速求解。

2 離散伴隨方程構(gòu)造的主要模塊

顯然，離散伴隨方程構(gòu)造核心是對(duì)Navier-Stokes方程右端殘差項(xiàng)進(jìn)行變分推導(dǎo)，涉及到對(duì)流項(xiàng)、人工黏性項(xiàng)、黏性項(xiàng)部分，以及邊界條件的處理。對(duì)各項(xiàng)的變分可以進(jìn)行手工推導(dǎo)，也可以借助自動(dòng)微分工具(如Tapenade、ADIFOR等)的后向模式來完成，前者的優(yōu)點(diǎn)是程序運(yùn)行效率較高、不依賴于第三方庫支持，缺點(diǎn)是工作繁瑣，容易出錯(cuò)；后者的特點(diǎn)是簡(jiǎn)捷方便，依賴于第三方庫支持、計(jì)算效率略低以及內(nèi)存需求偏大等問題。本文采用手工推導(dǎo)方式。

2.1 離散伴隨方程對(duì)流項(xiàng)的處理

離散伴隨方程對(duì)流項(xiàng)的構(gòu)造的主要依賴于空間離散格式、插值精度的選擇，不同的空間離散格式以及插值精度會(huì)產(chǎn)生不同的模板需求，尤其對(duì)于高精度格式來講，其無黏項(xiàng)的離散伴隨構(gòu)造將及其復(fù)雜。本文采用二階精度的中心格式，該格式構(gòu)造簡(jiǎn)單，在亞、跨、超聲速流場(chǎng)數(shù)值模擬中表現(xiàn)魯棒，在實(shí)際工程應(yīng)用較多。

對(duì)于l方向，式(8)的右端第二項(xiàng)中δR(j,k,l)可以表達(dá)為：

= (f1δξx)l+1/2+(ξxδf1)l+1/2+

(f2δξy)l+1/2+(ξyδf2)l+1/2+

(f3δξz)l+1/2+(ξzδf3)l+1/2-

[(f1δξx)l-1/2+(ξxδf1)l-1/2]-

[(f2δξy)l-1/2+(ξyδf2)l-1/2]-

其中，f1、f2、f3分別為無黏通量在笛卡爾坐標(biāo)系下三個(gè)方向的分量，對(duì)于定常問題，式(9)中忽略對(duì)幾何參數(shù)的變分后轉(zhuǎn)化為：

δR(j,k,l)=δfl+1/2-δfl-1/2

= (ξxδf1)l+1/2+(ξyδf2)l+1/2+(ξzδf3)l+1/2-

將上式代入式(9)并考慮以下關(guān)系式：

整理包含δWl的項(xiàng)，可以看到，對(duì)l單元伴隨無黏通量有貢獻(xiàn)的模板單元有l(wèi)-1，l，l+1，拓展到三維問題，對(duì)(j,k,l)單元伴隨無黏通量有貢獻(xiàn)的模板單元有：(j,k,l)，(j-1,k,l)，(j+1,k,l)，(j,k-1,l)，(j,k+1,l)，(j,k,l-1)，(j,k,l+1)共七個(gè)單元，如圖1所示。綜合式(8)、式(10)可以推導(dǎo)出離散伴隨方程的對(duì)流項(xiàng)：

圖1 離散伴隨無黏項(xiàng)模板單元Fig.1 Template unit of discrete adjoint inviscid term

2.2 離散伴隨方程人工黏性的處理

由于對(duì)單元(j,k,l)的伴隨人工黏性通量有貢獻(xiàn)的模板單元較多，因此，人工黏性的變分比較復(fù)雜，文中固定人工黏性系數(shù)前提下，給出了兩種處理方式，一種是直接變分推導(dǎo)，另外一種是對(duì)四階耗散項(xiàng)進(jìn)行簡(jiǎn)化處理，在梯度驗(yàn)證部分給出了兩種處理方式的可行性對(duì)比。

第一種方式～直接變分推導(dǎo)，首先僅考慮j方向：

δRD(j,k,l)=δfd,j+1/2-δfd,j-1/2

由式(12)可以看出，一維方向上對(duì)單元(j,k,l)的伴隨人工黏性通量有貢獻(xiàn)的模板單元有5個(gè)，推廣到三維問題，模板單元有13個(gè)，如圖2所示。

圖2 離散伴隨人工黏性項(xiàng)模板單元Fig.2 Template unit of discrete adjoint artificial viscous term

第二種方式～四階耗散項(xiàng)簡(jiǎn)化處理[18]，僅考慮j方向：

δRD(j,k,l)=δfd,j+1/2-δfd,j-1/2

上述簡(jiǎn)化，降低了推導(dǎo)難度，且虛網(wǎng)格單元的伴隨變量能夠利用邊界條件做到簡(jiǎn)單處理(見邊界條件部分)，由式(13)可以看出，一維方向上對(duì)單元(j,k,l)的伴隨人工黏性通量有貢獻(xiàn)的模板單元有3個(gè)，推廣到三維問題，模板單元有7個(gè)，如圖1所示。

2.3 離散伴隨方程黏性項(xiàng)的處理

離散伴隨方程黏性項(xiàng)的推導(dǎo)是最為復(fù)雜的一項(xiàng)，其核心是對(duì)速度導(dǎo)數(shù)項(xiàng)變分，如果直接采用完全NS方程進(jìn)行推導(dǎo)，將涉及速度導(dǎo)數(shù)交叉項(xiàng)，所需要的模板更多，推導(dǎo)將更加繁瑣，因此，本文采用黏性項(xiàng)薄層近似進(jìn)行變分。

曲線坐標(biāo)系下，采用如下薄層近似方式：

類似2.1部分的推導(dǎo)，從式(14)可以看出對(duì)單元(j,k,l)的伴隨黏性通量有貢獻(xiàn)的模板單元。

在式(14)的條件下，曲線坐標(biāo)系下的黏性通量可以表達(dá)為：

其中，

進(jìn)一步利用原始變量對(duì)守恒變量的轉(zhuǎn)換矩陣，求出對(duì)守恒變量的導(dǎo)數(shù)矩陣：

2.4 離散伴隨方程時(shí)間項(xiàng)隱式化

將上述推導(dǎo)結(jié)果進(jìn)行整合，并加入偽時(shí)間項(xiàng)可以得到離散伴隨主控方程：

Rc(λ)j,k,l-RD(λ)j,k,l-Rv(λ)j,k,l=0

對(duì)式(19)的迭代求解，可以采用顯式經(jīng)典四步龍格-庫塔推進(jìn)，也可以采用隱式時(shí)間推進(jìn)，這里我們將重點(diǎn)介紹LU-SGS方法，由于(19)在形式上與NS方程一致，因此，LU-SGS方法及其最大特征值分裂方法可以用于離散伴隨求解：

Qn+1=Qn+ΔQ

由于離散伴隨方程雅克比矩陣轉(zhuǎn)置的原因，上式中對(duì)應(yīng)的矩陣均需要進(jìn)行轉(zhuǎn)置處理，且無矩陣算法不再適用，右端項(xiàng)必須嚴(yán)格按照矩陣相乘進(jìn)行運(yùn)算，這是伴隨方程求解單步耗時(shí)、內(nèi)存需求高于NS方程的一個(gè)主要原因。流場(chǎng)時(shí)間推進(jìn)采用的隱式邊界條件在離散伴隨方程中依然可用：

3 求解離散伴隨方程的若干問題及處理方式研究與討論

對(duì)流項(xiàng)邊界條件、人工黏性項(xiàng)邊界條件、黏性項(xiàng)邊界條件、加速收斂技術(shù)、復(fù)雜項(xiàng)的簡(jiǎn)化方式、對(duì)離散伴隨方程求解來講至關(guān)重要，直接影響到梯度計(jì)算精度以及求解效率，為此，該部分將對(duì)上述關(guān)鍵環(huán)節(jié)進(jìn)行研究討論，探討對(duì)不同流動(dòng)問題的影響程度。

3.1 離散伴隨方程對(duì)流項(xiàng)的邊界條件

離散伴隨無黏通量計(jì)算在邊界處的模板單元與內(nèi)部點(diǎn)的個(gè)數(shù)不一樣，因此在邊界處的變分需要特殊處理，由于ξ,η,γ三個(gè)方向上的獨(dú)立性，對(duì)于ξ,η,γ=1,JKLDIM不同的邊界，可以分開處理，顯然，對(duì)于任意曲線坐標(biāo)方向，邊界單元無黏通量計(jì)算模板將減少一個(gè)，如圖3所示。

(a) J=1邊界

(b) J=JDIM邊界

以J=1邊界緊鄰的J=2單元通量計(jì)算為例，式(11)將改寫為(其他方向類似)：

式(25)中E,MBC分別對(duì)應(yīng)單位矩陣以及邊界條件矩陣，可以看出，離散伴隨無黏項(xiàng)的不同邊界條件變分，只需替換對(duì)應(yīng)邊界條件矩陣MBC，很容易實(shí)現(xiàn)不同邊界類型、以及內(nèi)外流伴隨之間的轉(zhuǎn)換、匹配，且這對(duì)于模塊化編程也十分有利。需要指出的是，由于離散伴隨無黏項(xiàng)的主導(dǎo)作用，該項(xiàng)邊界條件處理很大程度直接影響梯度的計(jì)算精度。

3.2 離散伴隨方程人工黏性的邊界條件

離散伴隨人工黏性通量計(jì)算模板單元較多，在邊界處的變分處理比較繁瑣，為此，而文中考察了2.2節(jié)的兩種處理方式，與2.1節(jié)一樣，其核心問題仍然是推導(dǎo)邊界條件矩陣MBC，由于推導(dǎo)過程較長(zhǎng)，本節(jié)僅給出了簡(jiǎn)化推導(dǎo)方式[19]，其他方式類似：

δRD(j,k,l)=δfd,j+1/2-δfd,j-1/2

圖4、圖5分別給出了某型飛機(jī)參數(shù)化及人工黏性處理方法對(duì)梯度計(jì)算精度影響比較 。

圖4 某型客機(jī)FFD參數(shù)化Fig.4 FFD lattice for wide body airplane

圖5 人工黏性處理方式對(duì)梯度計(jì)算精度影響比較Fig.5 Effect on gradient calculation of different artificial viscous processing method

3.3 離散伴隨方程黏性項(xiàng)的邊界條件

離散伴隨黏性通量的模板單元也將在邊界處發(fā)生變化。從式(16)的推導(dǎo)不難看出，黏性項(xiàng)的邊界條件需要考慮對(duì)虛網(wǎng)格的變分關(guān)系，實(shí)質(zhì)上仍然是推導(dǎo)邊界條件矩陣MBC。這樣，無論是無黏項(xiàng)、人工黏性項(xiàng)還是物理黏性項(xiàng)的邊界條件均轉(zhuǎn)化為一個(gè)矩陣推導(dǎo)，大大簡(jiǎn)化了程序設(shè)計(jì)框架，式(26)～式(29)給出了幾類典型的邊界條件矩陣。

廣義對(duì)稱邊界與無黏物面條件矩陣：

黏性物面條件矩陣：

超聲速入流/出流邊界條件矩陣：

其他類型邊界條件矩陣同樣可以依據(jù)邊界類型進(jìn)行求導(dǎo)。

3.4 離散伴隨方程的時(shí)間推進(jìn)

不同時(shí)間推進(jìn)方式，同樣也對(duì)離散伴隨方程產(chǎn)生不同的收斂效果，比如GMRES、BLU-SGS時(shí)間推進(jìn)方式，綜合考慮內(nèi)存需求與計(jì)算效率，文中進(jìn)行了BLU-SGS與LU-SGS收斂歷程以及計(jì)算結(jié)果的對(duì)比，在伴隨方程求解中依然有效。

圖6、圖7分別給出了不同的時(shí)間推進(jìn)方式對(duì)收斂速度及計(jì)算結(jié)果的對(duì)比。

圖6 LU-SGS、BLU-SGS收斂速度比較Fig.6 Convergence rate comparison of LU-SGS and BLU-SGS

圖7 LU-SGS、BLU-SGS第一伴隨變量計(jì)算結(jié)果Fig.7 The calculation result of fist adjoint variable of LU-SGS and BLU-SGS

3.5 離散伴隨方程的并行化求解

與流場(chǎng)并行計(jì)算一樣，離散伴隨方程求解時(shí)，并行機(jī)制依然采用單元數(shù)衡量的負(fù)載平衡、對(duì)等式計(jì)算以及MPI消息傳遞模式，對(duì)于伴隨方程來講，依賴于求解器的構(gòu)架，通過MPI進(jìn)行傳遞的信息可以是雅克比矩陣，也可以是伴隨變量本身。本文求解器采用了多塊對(duì)接網(wǎng)格技術(shù)，與對(duì)接面邊界信息一樣，MPI傳遞的信息是各個(gè)進(jìn)程分割面上的兩層虛網(wǎng)格上的伴隨變量信息，這樣離散伴隨求解的并行效率特性與流場(chǎng)的基本一致，由于存在矩陣運(yùn)算，略低于流場(chǎng)并行效率。圖8給出了NS方程與離散伴隨方程并行效率比較。

圖8 離散伴隨/NS方程并行效率Fig.8 Parallel efficiency of discrete adjoint and NS equation

3.6 離散伴隨方程求解的加速收斂技術(shù)

在流場(chǎng)計(jì)算中采用的當(dāng)?shù)貢r(shí)間步長(zhǎng)、多重網(wǎng)格、網(wǎng)格序列法等加速收斂技術(shù)，在離散伴隨方程求解中使用仍然能夠起到加速作用。無論是多重網(wǎng)格方法，還是網(wǎng)格序列法，均要涉及粗網(wǎng)格殘差計(jì)算方法，本文在粗網(wǎng)格上的殘差計(jì)算采用二階格式且不考慮黏性項(xiàng)，而對(duì)于伴隨變量本身的插值、限制方式仍與流場(chǎng)變量一致；當(dāng)?shù)貢r(shí)間步長(zhǎng)取值與流場(chǎng)計(jì)算相同，不同的是，求解離散伴隨方程時(shí)是基于流場(chǎng)收斂解進(jìn)行的，因此，穩(wěn)定性更好，CFL可以取得更大(圖9)。圖10為超聲速低聲爆客機(jī)第一伴隨變量云圖，在黏性流場(chǎng)計(jì)算時(shí)CFL≤10.0取值基本保證穩(wěn)定收斂，離散伴隨計(jì)算時(shí)CFL=50依然能夠穩(wěn)定收斂。

圖9 網(wǎng)格序列法加速收斂(CFL=50 無黏伴隨)Fig.9 Convergence accelerated method based on grid sequence method(Inviscid adjoint)

圖10 超聲速低聲爆客機(jī)第一伴隨變量云圖 (目標(biāo)函數(shù)：阻力CFL=50)Fig.10 The first adjoint variable contour of supersonic low sonic boom airplane

4 典型工程算例測(cè)試

求解方法對(duì)梯度計(jì)算的影響經(jīng)過討論之后，進(jìn)一步采用本文的方法開展不同流動(dòng)問題梯度計(jì)算精度驗(yàn)證，采用了兩種構(gòu)型：某型客機(jī)全機(jī)巡航構(gòu)型和外壓式超聲速進(jìn)氣道。

4.1 外部黏性繞流問題——某型客機(jī)全機(jī)巡航構(gòu)型

以某型某型客機(jī)全機(jī)巡航構(gòu)型外部繞流為算例，主要部件包含機(jī)翼、機(jī)身、掛架、短艙內(nèi)外涵道、平尾以及立尾。半模網(wǎng)格劃分為526塊，網(wǎng)格規(guī)模2500萬量級(jí)，如圖11所示。采用SST湍流模型，128核進(jìn)行并行計(jì)算。

圖12給出了黏性離散伴隨方程的收斂歷程，可以看出文中的求解方法對(duì)于復(fù)雜外形也非常穩(wěn)定。圖13、圖14分別為物面第一伴隨變量云圖與參數(shù)化示意圖。圖15給出了任意選取的幾個(gè)控制頂點(diǎn)的導(dǎo)數(shù)值對(duì)比，梯度幅值以及梯度方向一致，平均誤差為6%，可以看出，在外部黏性擾流問題中，基于薄層近似的黏性離散伴隨導(dǎo)數(shù)計(jì)算精度較高，完全滿足工程氣動(dòng)設(shè)計(jì)要求。

圖11 表面網(wǎng)格分布Fig.11 Grid distribution on surface

圖12 離散伴隨方程收斂歷程Fig.12 Convergence process of discrete adjoint equation

圖13 物面第一伴隨變量云圖Fig.13 Contour of the first adjoint variable on surface

圖14 機(jī)翼參數(shù)化FFD lattice示意圖Fig.14 FFD lattice for wing parameterized

圖15 黏性離散伴隨梯度計(jì)算與差分對(duì)比Fig.15 Comparison of gradient calculation of viscous discrete adjoint and difference

4.2 內(nèi)外流一體化問題——外壓式無附面層隔道超聲速進(jìn)氣道(DSI進(jìn)氣道)

該算例為超聲速無附面層隔道進(jìn)氣道內(nèi)外流一體化數(shù)值模擬，涉及超聲速入口、超聲速出口、黏性物面以及風(fēng)扇入口(質(zhì)量流流場(chǎng)出口)四類邊界條件。

需要指出的是，流場(chǎng)計(jì)算直接采用質(zhì)量出口邊界條件，帶來收斂速度慢、邊界條件矩陣推導(dǎo)困難等問題，因此本文采用了壓力特征出口—質(zhì)量流調(diào)節(jié)的方法，即給定初始反壓，計(jì)算程序依據(jù)流量特征進(jìn)行自動(dòng)調(diào)節(jié)，該方法可以提高管道內(nèi)“激波外推”速度，加速收斂。流場(chǎng)收斂后風(fēng)扇入口流場(chǎng)仍然保持亞聲速壓力特征出口邊界，有利于其邊界條件矩陣推導(dǎo)。

圖16給出了設(shè)計(jì)狀態(tài)進(jìn)氣道波系狀態(tài)分布，該模型的初始鼓包依據(jù)三維乘波理論設(shè)計(jì)，圖17給出了鼓包FFD參數(shù)化示意圖，沿進(jìn)氣方向選取5個(gè)控制點(diǎn)，進(jìn)行總壓恢復(fù)系數(shù)梯度測(cè)試，這里給出總壓對(duì)守恒變量的變分推導(dǎo)方法：

圖18給出了鼓包上任意選取幾個(gè)控制點(diǎn)，計(jì)算風(fēng)扇入口總壓對(duì)控制點(diǎn)的導(dǎo)數(shù)，由于內(nèi)外流一體化問題流動(dòng)黏性效應(yīng)較強(qiáng)，而文中黏性離散伴隨計(jì)算采用了薄層近似，本文將對(duì)此進(jìn)一步開展研究，消除這些因素帶來的誤差。綜上所述，文中采用的方法誤差較為明顯，但依然較好地反映了梯度幅值變化趨勢(shì)，且梯度方向一致，能夠應(yīng)用于氣動(dòng)優(yōu)化。

圖16 無附面層隔道進(jìn)氣道波系分布Fig.16 Wave system distribution of DSI

圖17 DSI進(jìn)氣道鼓包參數(shù)化Fig.17 The parameterized DSI inlet bump

圖18 黏性離散伴隨梯度計(jì)算與差分對(duì)比Fig.18 Comparison of gradient calculation of viscous discrete adjoint and difference method

5 結(jié) 論

基于自主研發(fā)的大規(guī)模并行結(jié)構(gòu)化網(wǎng)格RANS求解器PMB3D，開展了黏性離散伴隨方程構(gòu)造、求解方法的研究與討論：

1) 以矩陣的形式體現(xiàn)離散伴隨方程的邊界條件，有利于程序的簡(jiǎn)捷化、模塊化。

2) 人工黏性的簡(jiǎn)化處理雖然降低了的梯度求解精度，但仍然能夠較好地反映梯度變化趨勢(shì)，可以應(yīng)用于氣動(dòng)優(yōu)化。

3) 結(jié)構(gòu)化求解器網(wǎng)格中，無黏項(xiàng)、人工黏性項(xiàng)的離散推導(dǎo)具有方向獨(dú)立性；薄層假設(shè)條件下，黏性項(xiàng)的離散推導(dǎo)能夠大大簡(jiǎn)化推導(dǎo)過程以及實(shí)現(xiàn)程序模塊化。

4) 由于存在大量的矩陣運(yùn)算，離散伴隨方程求解計(jì)算效率、并行效率均低于NS方程求解。

5) 薄層假設(shè)條件在外流問題梯度計(jì)算中精度較高，在內(nèi)流問題計(jì)算中梯度精度降低，但可以較好地反映梯度變化趨勢(shì)。

致謝:在本文的研究工作中，肖中云、牟斌、李彬、唐靜、賈洪印等同志在大規(guī)模并行塊隱式時(shí)間推進(jìn)方法、進(jìn)氣道計(jì)算等方面給予了有益的建議，在此表示感謝！

[1]Jameson A. Aerodynamic design via control theory[J]. Journal of Scientific Computing, 1988, 3: 233-260

[2]Giles M B, Duta M C. Algorithm developments for discrete adjoint methods[J]. AIAA Journal, 2003, 41(2): 198-205

[3]Carpentieri G. An adjoint-based shape-optimization method for aerodynamic design[D]. Delft: Technische Universiteit, 2009

[4]Amoignon O, Berggren M. Adjoint of a median-dual finite-volume scheme application to transonic aerodynamic shape optimization[R]. Technical Report 2006-3, Uppsala University, 2006

[5]Reuther J. Aerodynamic shape optimization using control theory[D]. Ph.D. Dissertation, University of California, Davis, Davis, CA, June 1996

[6]Thomas A Zang, et al. Multidisplinary design optimization techniques: implications and opportunities for fluid dynamics research[C]//30thAIAA Fluid Dynamics conference, Norfolk, VA. June 28-July, 1999

[7]Nielsen E J, Anderson W K. Recent improvements in aerodynamic design optimization on unstructrued meshes[J]. AIAA Journal, 2002, 40(6): 1155-1163

[8]Dwight R P, Brezillon J. Effect of various approximations of the discrete adjoint on gradient-based optimiza-tion[R]. AIAA-2006-0690, 2006

[9]Carrier G, Destarag D, et al. Gradient-based aerodynamic optimization with the elsA software[R]. AIAA 2014-0568

[10]Qin N, Wong W S, Moigne A L. Three-dimensional contour bumps for transonic wing drag reduction[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2008, 222(5): 619-629

[11]Zuo Yingtao, Gao Zhenghong, Zhan Hao. Aerodynamic design method based on N-S equations and discrete adjoint approach[J]. Acta Aerodynamica Sinica, 2009, 27(1): 67-72. (in Chinese)左英桃, 高正紅, 詹浩. 基于N-S方程和離散共軛方法的氣動(dòng)設(shè)計(jì)方法研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2009, 27(1): 67-72

[12]Xiong Juntao, Qiao Zhide, Yang Xudong, et al. Optimum aerodynamic design of transonic wing based on viscous adjoint method[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(2): 281-285.(in Chinese)熊俊濤, 喬志德, 楊旭東, 等. 基于黏性伴隨方法的跨聲速機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2007, 28(2): 281-285

[13]Qu Kun, Li Jichao, Cai Jinsheng. Method of linearizing computational fluid dynamics model and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3218-3227.(in Chinese)屈崑, 李記超, 蔡晉生. CFD數(shù)學(xué)模型的線性化方法及其應(yīng)用[J]. 航空學(xué)報(bào), 2015, 36(10): 3218-3227

[14]Zhang Chaolei, Li Haitao, Feng Zhenping. Aerodynamic optimization design of turbomachinery cascade based on discrete adjoint method[J]. Journal of Engineering Themophics, 2012, 33(1): 47-50.(in Chinese)張朝磊, 等. 基于離散伴隨方法的透平葉柵氣動(dòng)優(yōu)化[J]. 工程熱物理學(xué)報(bào), 2012, (33): 47-50

[15]Gao Yisheng, Wu Yizhao, Xia Jian. A discrete adjoint-based approach for airfoil optimization on unstructured meshes[J]. Acta Aerodynamica Sinica, 2013, 31(2): 244-249.(in Chinese)高宜勝, 伍貽兆, 夏健. 基于非結(jié)構(gòu)網(wǎng)格離散型伴隨方法的翼型優(yōu)化[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(2): 244-249

[16]Li Bin, Deng Youqi, Tang Jing, et al. Discrete adjoint optimization method for 3d unstructured grid[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 674-686. (in Chinese)李彬, 鄧有奇, 唐靜, 等. 基于三維非結(jié)構(gòu)混合網(wǎng)格的離散伴隨優(yōu)化方法[J]. 航空學(xué)報(bào), 2014, 35(3): 674-686

[17]Liu Gang, Xiao Zhongyun, Wang Jiantao, Liu Fan. Numerical simulation of missile air-launching process under rail slideway constraints[J]. Acta Aerodynamica Sinica, 2015, 33(02): 192-197.(in Chinese)劉剛, 肖中云, 王建濤, 等. 考慮約束的機(jī)載導(dǎo)彈導(dǎo)軌發(fā)射數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2015, 33(02): 192-197

[18]Nadarajah Siva K, Jameson Antony. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization[C]//38thAIAA Aerospace Sciences Meeting and Exhibit Reno-NV-Jan. 10-13-2000

[19]Chen R F, Wang Z J. Fast, block lower-upper symmetric gauss-seidel scheme for a artributrary grids[J]. AIAA Journal, 2000, 38(12): 2238-2245.

Investigation of gradient computation based on discrete adjoint method

HUANG Jiangtao1, LIU Gang1, ZHOU Zhu1,*, GAO Zhenghong2, HUANG Yong1
(1. Computational Aerodynamics Institute of China Aerodynamics Research and Development Center, Mianyang 621000, China; 2. Northwestern Polytechnical University, Xi’an 710072, China)

Viscous discrete adjoint equations and the corresponding solving method are studied and discussed based on PMB3D, i.e., a parallelized in-house CFD code for multi-block structured grid. Firstly, the discrete adjoint gradient strategy is briefly introduced. Secondly, we fully describe the variation derivation and the contribution of inviscid, artificial viscosity, and viscous parts to the discrete adjoint equations. Thirdly, the implementation of the inviscid part and the viscous boundary conditions are studied, and the simplification methods for variation derivation are discussed. Two typical simulations are respectively carried out on a wide-body aircraft and an external compression supersonic inlet to analyze the influences of these simplification methods on the solution precision. Finally, parallel solution, time integration, and convergence acceleration are discussed and validated. Numerical simulation demonstrates that the present discrete adjoint gradient strategy and equations contribute to a convenient programming and modularization. Moreover, the precision of the gradient solution is qualified for aerodynamic optimization design. The present study can be recognized as a useful reference for further study on discrete adjoint method.

discrete adjoint; aerodynamic optimization; gradient information; boundary condition; parallelized computation

0258-1825(2017)04-0554-09

2017-04-02;

2017-06-28

國家自然科學(xué)基金(11402288)；國家重點(diǎn)研發(fā)計(jì)劃“數(shù)值飛行器原型系統(tǒng)”(2016YFB0200704)

黃江濤(1982-)，男，副研究員，研究方向：飛行器氣動(dòng)外形多學(xué)科優(yōu)化與計(jì)算空氣動(dòng)力學(xué). E-mail: hjtcyf@163.com

周鑄(1973-)，男，重慶人，研究員，研究方向：飛行器氣動(dòng)設(shè)計(jì)與計(jì)算空氣動(dòng)力學(xué). E-mail: zhouzhu@tom.com

黃江濤, 劉剛, 周鑄, 等. 基于離散伴隨方程求解梯度信息的若干問題研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(4): 554-562.

10.7638/kqdlxxb-2017.0064 HUANG J T, LIU G, ZHOU Z, et al. Investigation of gradient computation based on discrete adjoint method[J]. Acta Aerodynamica Sinica, 2017, 35(4): 554-562.

V211.3

A doi: 10.7638/kqdlxxb-2017.0064