黃江濤，周鑄,*，劉剛，高正紅，黃勇，王運濤

1.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000 2.西北工業(yè)大學 翼型葉柵空氣動力學國防科技重點實驗室， 西安 710072

基于伴隨方程的氣動優(yōu)化以其獨有的優(yōu)勢，在氣動設計等領域扮演著重要角色，也是國內外空氣動力學研究機構一個重要的研究方向，而基于變分思想的氣動/結構多學科伴隨優(yōu)化方法也開始在工程領域發(fā)揮重要作用，例如考慮氣動彈性變形的柔性機翼若采用基于差分的梯度優(yōu)化以及進化算法開展多學科多目標優(yōu)化，其計算量非常龐大，甚至難以忍受，設計效率極為低下。另一方面，未來飛機發(fā)展的一個重要方向是重量較輕的復合材料結構柔性機翼設計(如B787、B747-8等寬體客機)，此時氣動、結構耦合效果將更加明顯，基于多學科耦合伴隨系統(tǒng)的靈敏度分析在綜合設計上具有更加突出的優(yōu)勢。針對該研究方向，國外在多學科耦合伴隨[1]方面開展了系列的研究，例如密歇根大學Martins教授的多學科優(yōu)化設計(MDO)團隊，基于結構網格CFD求解器以及有限元方法發(fā)展了氣動/結構延遲耦合伴隨(Lagged Coupled Adjoint,LCA)方法，實現(xiàn)了氣動結構一體化設計[2-3]，斯坦福大學Jameson團隊基于耦合伴隨系統(tǒng)進行了機翼平面形狀與剖面的氣動結構多學科優(yōu)化[4]，德宇航Abu-Zurayk和Brezillon基于非結構化求解器TAU[5]，以及法宇航Marcelet等基于CFD代碼elsA也發(fā)展了LCA優(yōu)化方法[6-7]。

國內在流場單學科離散伴隨方程求解器以及基于代理模型的氣動/結構多學科優(yōu)化方面取得了一定的進展[8-14]。大多工作局限于單學科伴隨系統(tǒng)以及基于進化算法的優(yōu)化，在多學科耦合伴隨系統(tǒng)自主研發(fā)方面研究基礎比較薄弱。由于氣動/結構耦合伴隨系統(tǒng)具有計算代價小，梯度計算量與各個學科設計變量個數(shù)均無關等優(yōu)點，且通過耦合伴隨方程的求解能夠快速計算出氣動、結構等學科關心的各個目標函數(shù)對公共設計變量以及獨立設計變量的導數(shù)，將在未來多學科優(yōu)化領域發(fā)揮重要作用。

由于耦合伴隨系統(tǒng)中交叉學科導數(shù)項的具體推導方法，耦合方程組對應左端的各個大型稀疏矩陣的求解方式和變分簡化處理方式，以及學科之間物理場信息、伴隨變量交換方式直接影響了多學科變分的簡捷性、耦合系統(tǒng)計算效率與梯度信息的計算精度。因此，本文將對以上幾個關鍵環(huán)節(jié)的變分推導、各個學科伴隨方程的求解方法進行詳細系統(tǒng)的研究與改進?；谡n題組自主研發(fā)的大型并行結構化網格雷諾平均Navier-Stokes (RANS)求解器PMB3D (Paralleized Muti-Block solver for 3-Dimensional flow)[15]以及流固耦合代碼FSC3D(Fluid-Solid Coupled platform for 3-Dimensional case)首先建立了飛行器靜氣動彈性數(shù)值模擬技術，進一步基于并行化伴隨方程求解器PADJ3D (Paralleized ADJoint solver for 3-Dimensional flow)[16]，開展多學科耦合伴隨系統(tǒng)的構建與求解研究，希望能夠為多學科耦合伴隨系統(tǒng)的研究人員提供有價值的參考。

1 耦合系統(tǒng)中的學科分析代碼

1.1 并行化CFD代碼PMB3D

PMB3D是課題組自主研發(fā)的大型并行CFD代碼，求解任意曲線坐標系下的Navier-Stokes方程實現(xiàn)流場精細化數(shù)值模擬,即

(1)

PMB3D具備S-A一方程、剪切應力輸運(SST)兩方程湍流模型以及Langtry-Menter轉捩預測模型， Roe、Vanleer等多種空間離散方法，MUSCL(Monotonic Upwind Scheme for Conservation Laws)迎風插值與WENO(Weighted Essentially Non-Oscillatory )格式，隱式時間推進，多重網格技術，支持多塊對接、拼接以及重疊網格技術，以及基于MPI(Message Passing Interface)通信協(xié)議的大規(guī)模并行計算能力，廣泛應用于常規(guī)氣動力計算、多體分離、空中加油安全性分析、進排氣系統(tǒng)、旋轉部件氣動特性計算以及火箭發(fā)射、級間分離等領域[15]。

1.2 流固耦合代碼FSC3D

FSC3D代碼采用虛功原理與徑向基函數(shù)(RBF)插值技術[17]實現(xiàn)CFD邊界面元氣動力向結構點、結構點位移向CFD物面格點高精度插值，基于RBF方法的插值技術數(shù)學模型表示為

(2)

式中：αi為插值系數(shù);φ(·)為基函數(shù)。

本文對p(x)的選取采用多項式方法，即

式中：αi(i=0,1,…,n)為插值系數(shù)；xi,n為中心點坐標;M為徑向基中心矢量維數(shù)。

1.3 靜氣動彈性數(shù)值模擬代碼

采用LDLT方法[18]求解結構靜力學方程：

Kd=F

(3)

式中：K為剛度矩陣；d為結構節(jié)點在6個自由度的變形位移；F為作用在節(jié)點上的氣動力矩陣。在完成位移耦合即氣動表面網格更新后，本文采用徑向基結合無限插值技術進行空間網格變形，由上述RBF技術依據(jù)物面網格頂點的變形量構建RBF精確插值模型, 操作空間網格Block頂點進行變形插值, 保證網格整體空間拓撲一致性，進一步利用TFI(TransFinite Interpolation)無限插值進行Block中ξ、η、γ3個方向的Edge、Face以及內點更新[19],即

(4)

式中：dx和s、t、u分別為位移以及3個方向上的邏輯坐標；NI、NK、NJ分別代表3個方向的網格維度;Sξ,η,γ、tξ,η,γ和uξ,η,γ代表3個曲線坐標方向。

為了加速靜氣動彈性數(shù)值模擬的收斂速度，本文采用流固耦合計算與CFD求解器共用模塊的方式，即在流場迭代至指定步數(shù)，利用不完全收斂解進行彈性變形計算，直至流場收斂，如圖1所示。

圖1 共用solver模塊方法與傳統(tǒng)氣動彈性計算方法流程對比Fig.1 Comparison of process of common solver module method with traditional aeroelasticity calculation method

2 流場伴隨方程簡介

對于氣動優(yōu)化設計的最小化問題：

minI(w,X,D)

(5)

式中：I可以是升、阻力、力矩、流量、壓比等參數(shù)；w、X、D分別為流場守恒變量、網格變量、設計變量;在流場收斂解的條件下，殘差約束R(w,X,D)=0，引入拉格朗日算子λ構造目標函數(shù)，即

L=I+λTR

(6)

對式(6)進行求導，可得

(7)

(8)

式(8)就是流場伴隨方程，通過隱式迭代方法求解拉格朗日算子λ之后，可以通過式(9)進行目標函數(shù)梯度信息快速求解,即

(9)

(10)

本文采用二階精度的中心格式、人工黏性以及SST兩方程湍流模型，手工推導構造Navier-Stokes方程離散伴隨方程，并進一步引入虛擬時間項：

(11)

式中：Vj,k,l為網格體積；Rc(λ)、RD(λ)、Rv(λ)j,k,l分別代表伴隨方程的對流項、人工黏性項以及物理黏性項。對式(11)的迭代求解，采用LU-SGS方法隱式時間推進，離散伴隨方程各項的空間離散采用矩陣形式的物理邊界條件，凍結湍流黏性系數(shù)，時間推進采用的隱式邊界條件為

Δλ*=0, Δλ=0

(12)

式中:λ*為隱式推進的中間項。

黏性項采用薄層近似，曲線坐標系下的黏性通量可以表達為

(13)

式中:各項表達式具體含義可參考文獻[16]。

與流場并行計算一樣，離散伴隨方程求解時，并行機制依然采用單元數(shù)衡量的負載平衡、對等式計算以及MPI消息傳遞模式，對于伴隨方程來講，依賴于求解器的構架，通過MPI進行傳遞的信息可以是雅克比矩陣，也可以是伴隨變量本身。本文求解器采用了多塊對接網格技術，與流場對接面邊界信息一樣，MPI傳遞的信息是各個進程中分割面上的兩層虛網格上的伴隨變量信息。

3 考慮結構靜力學方程的耦合伴隨方程

3.1 耦合伴隨方程

參照第2節(jié)流場變分思想，綜合考慮最小化目標函數(shù)、流場收斂解約束以及氣動彈性平衡約束，有

minI(w,X,D)

(14)

式中：Ra、Rs分別為流場殘差與結構靜力學殘差；K、d、F分別為結構剛度矩陣、變形位移以及載荷分布。對式(14)引入拉格朗日算子ψ可以構造目標函數(shù)：

(15)

(16)

令式(16)第1項為0，獲取氣動/結構多學科耦合伴隨方程為

(17)

即

(18)

3.2 延遲耦合伴隨方程

(19)

式(18)被稱之為延遲耦合伴隨方程，可以看出，通過延遲伴隨變量的引入，式(19)已經解耦，各個學科之間的影響通過方程右端的強迫項來實現(xiàn)，不同學科方程之間可以進行松耦合迭代，大幅度降低了方程求解難度。

其中，流場伴隨方程強迫項為

(20)

(21)

式中:wρ、wρ u、wρ v、wρ E和wρ w對應守恒變量的展開表達式。

由結構靜力學方程(結構節(jié)點個數(shù)為Ns),有

(22)

上述矩陣維度取決于結構節(jié)點的維度以及CFD格心單元的維度，ψs的維度取決于有限元模型的網格類型，對于板殼單元，ψs為6Ns×1維度的矢量，?Rs/?w的轉置矩陣As為5Na×6Ns矩陣。因此,結構殘差對流場變量的導數(shù)?Rs/?w通過有限差分求出，解析表達式推導方式依賴于所采用的流固耦合方法?；诒疚牡奶摴υ砹鞴恬詈戏绞剑旅娼o出具體的推導。

由于FSC3D代碼采用虛功原理實現(xiàn)氣動力向結構載荷插值[18]，x方向(其他方向表達形式完全一致)結構載荷與氣動力插值可以表達為矩陣與向量形式,即

(23)

(24)

式(24)直接對守恒變量求導比較困難，可以變換為

(25)

代入耦合伴隨方程，可以獲得流場伴隨方程強迫項最終表達形式。

同理，結構伴隨方程強迫項為

(26)

上述矩陣維度同樣取決于結構節(jié)點的維度以及CFD格心單元的維度，ψa為5Na×1維度的矢量，?Ra/?d的轉置矩陣Aa為6Ns×5Na矩陣。流場殘差對結構位移的導數(shù)?Ra/?d求解方式同樣依賴于流固耦合方法以及變形網格方法，本文采用RBF_TFI變形網格方法，即

(27)

式中：X、Xsurf分別為空間網格與物面網格。從式(27)右端項中可以看出，流場殘差對結構位移的導數(shù)?Ra/?d的解析表達形式依賴于所采用的動網格方法以及流固耦合技術，即第1項可以利用自動微分進行求導，對于彈簧法、徑向基類型的動網格來講，第2項可以推導出變形網格矩陣KD，第3項取決于結構位移向氣動表面網格的插值方式，由于文中采用了RBF插值技術，很容易手工推導出位移向表面網格的變換矩陣KR。需要指出的是，為提高結構化變形網格的運算效率，文中采用了RBF_TFI混合算法，此時很難顯式地推導出矩陣KD，該項的求導在文中由簡單的單側有限差分代替，由于RBF_TFI混合算法具有極高的效率，差分計算量可以忽略。

延遲伴隨方程第1項強迫項為大型矩陣與矢量相乘，將其具體表達式展開，下標表達式為

結構伴隨方程強迫項的下標表達式為

(28)

綜合強迫項的整理，可以推導出耦合延遲伴隨方程偽時間-殘差表達形式為

(29)

式(29)流場伴隨方程與結構伴隨方程一步可以利用式(11)中的LU-SGS隱式迭代或雅克比迭代方法進行求解。

觀察式(29)結構伴隨方程，其本質上是大型稀疏矩陣求解，對于該類矩陣求解，常用的LU分解以及高斯消去法帶來代價較大、在針對較多自由度結構有限元大型矩陣求解計算效率偏低等問題。綜合考慮計算效率與存儲量，Cholesky 分解算法是一個較為合理的選擇，Cholesky 分解算法中主要包含LLT以及LDLT方法[18]，前者要求矩陣必須正定，而后者沒有要求，魯棒性更強，相對于迭代方法以及傳統(tǒng)求逆方法，LDLT的計算效率較高，無需迭代且求解結果更為精確，因此本文采用并行LDLT方法代替迭代方法進行結構伴隨方程求解，對結構伴隨方程矩陣項進行LDLT分解，即

分解出Lij和dii元素后，結構伴隨方程求解可以通過三步掃描快速完成，即

(30)

式中：?和φ均為中間變量。

綜合以上求解過程，可以看出延遲伴隨方程的解耦求解方式，能夠充分利用原有的求解體系，并不破壞原有程序的基本框架，由此可以直觀地獲取LCA系統(tǒng)各個分析模塊的組裝關系與工作流程，如圖2所示。

對于延遲方程組式(19)中結構伴隨方程的右端項，其含義是計算目標函數(shù)對結構位移的直接導數(shù)以及流場殘差對結構位移導數(shù)與流場延遲伴隨變量的乘積，觀察單一流場學科目標函數(shù)導數(shù)的計算公式：

圖2 LCA方法中各個學科分析模塊組裝與流程Fig. 2 Disciplines assembly and work flow of LCA method

(31)

不難發(fā)現(xiàn),兩者在表達形式上完全一致，僅僅在自變量擾動上不同，前者是結構位移，后者是參數(shù)化控制頂點，因此，結構伴隨方程右端項計算與導數(shù)計算可以由同一模塊來完成，只需將子程序的實參由參數(shù)化控制頂點P_ffd替換為結構位移CSD_displacement即可，如圖3和圖4所示。

圖3 流場單學科目標函數(shù)導數(shù)計算模塊Fig.3 Objective function derivative module of flow field single discipline

圖4 結構伴隨方程右端項計算模塊Fig.4 Computational module for right end item of structural adjoint equation

進一步利用鏈式求導展開，即

(32)

可以看出，?X/?Xsurf的推導方式取決于所采用的變形網格方法，?Xsurf/?d的推導方式取決于所采用的流固耦合方法。

3.3 延遲耦合伴隨方程求解與梯度計算

對式(18)加入偽時間項進行迭代求解，流場延遲伴隨方程求解采用LU-SGS時間推進，結構延遲伴隨方程求解采用LDLT方法進行求解，間隔指定的子迭代步數(shù)進行一次延遲伴隨變量傳遞以及不同學科殘差導數(shù)矩陣運算，實現(xiàn)強迫項的交換。

流場伴隨方程求解采用單元數(shù)衡量的負載平衡，對等式計算以及MPI消息傳遞模式，本文求解器采用了多塊對接網格技術，所以MPI傳遞的信息是各個進程分割面上的兩層虛網格上的伴隨變量。求解過程中流場延遲伴隨方程的邊界條件仍然采用矩陣形式[16]，結構伴隨方程的相關導數(shù)矩陣由流固耦合程序FSC3D變分計算。求出延遲伴隨變量ψ=[ψaψs]T，代入導數(shù)求解公式，可以獲取目標函數(shù)的梯度為

(33)

式(33)的矩陣表達形式為

對各項進行進一步展開，有

式中：Swall、G分別代表物面網格以及空間網格。不難看出?G/?Swall的表達形式取決于所采用的變形網格方法，?Swall/?X的推導方式取決于所采用的參數(shù)化建模方法。

4 學科分析耦合系統(tǒng)的驗證

采用HIRENASD(HIgh REynolds Number Aero-Structural Dynamics)模型[20]進行靜氣動彈性計算，并與試驗數(shù)據(jù)進行對比，驗證耦合分析手段的正確性。

圖5給出了網格分布以及典型計算狀態(tài)(Ma=0.8，Re=7.0×106)的彈性變形示意圖；圖6給出了剛性與氣動彈性變形的機翼展向66%、95%站位壓力分布與風洞試驗的對比，其中,Cp為壓力系數(shù),x/c為展弦比。可以看出彈性變形數(shù)值模擬結果與試驗數(shù)據(jù)更為吻合，驗證了氣動彈性手段的可靠性，可以為多學科耦合伴隨系統(tǒng)數(shù)值模擬提供準確的多物理場信息。

圖5 CFD網格分布與彈性變形Fig.5 CFD grid distribution and aeroelastic deformation

圖6 壓力分布與風洞試驗對比Fig.6 Pressure distribution compared with wind tunnel test

5 典型工程算例測試

基于CRM (Common Research Model)標準算例構型，進行本文的方法梯度計算精度以及收斂特性驗證。

主要部件包含機翼、機身、平尾以及立尾。網格劃分為290塊，半模網格規(guī)模為1 260萬，如圖7所示。采用中心格式、SST湍流模型、LU-SGS時間推進以及多重網格加速收斂技術，64核并行計算。

圖8給出了基于非均勻有理B樣條(NURBS)基函數(shù)的自由式變形 (FFD)參數(shù)化示意圖，共采用200個控制頂點實現(xiàn)機翼氣動外形參數(shù)化建模；為方便耦合系統(tǒng)的測試，結構有限元建模采用簡化板殼單元，主要結構單元包含梁、翼肋等，結構參數(shù)化方法采用FFD技術，實現(xiàn)梁、翼肋的高度和寬度變化，有限元模型節(jié)點數(shù)為490個，如圖9和圖10所示。圖11給出了流固耦合與流場求解模塊共用一個模塊的殘差以及固定升力系數(shù)下氣動力的收斂歷程，每500步進行一次流固耦合模塊調用，每10步進行一次迎角調整，經過8次氣動彈性靜力學方程求解，5 000步左右的迭代，氣動力基本趨于收斂，與單一流場計算相比，計算量增加很小。

圖12給出了典型耦合步的彈性變形及其局部視圖。圖13給出了多學科耦合伴隨系統(tǒng)殘差收斂歷程，可以看出耦合系統(tǒng)的殘差收斂穩(wěn)定，且耦合產生的跳躍殘差也隨著迭代過程下降，這說明隨著迭代的收斂，多學科延遲耦合的解與原始耦合伴隨系統(tǒng)趨于一致；圖14給出了結構伴隨變量殘差隨迭代步數(shù)的收斂歷程，可以看出經過6次的強迫項耦合計算，結構伴隨方程的解基本收斂。進一步對耦合伴隨梯度的計算精度與差分方法進行對比，差分步長為0.000 1。

圖7 CFD表面網格分布Fig.7 CFD gird distribution on wall

圖8 自由式變形參數(shù)化Fig.8 Free-form deformation lattice for parameterization

圖9 蒙皮結構有限元網格Fig.9 Finite element grid of skin structure

圖10 梁、肋等有限元網格Fig.10 Finite element grid of beam and rib

圖11 氣動彈性收斂歷程Fig.11 Process of aeroelastic convergence

圖12 典型耦合步長彈性變形及局部視圖Fig.12 Elastic deformation with typical coupling step and local view of elastic deformation

圖13 耦合伴隨系統(tǒng)殘差收斂歷程Fig.13 Process of residual convergence of coupled adjoint system

求解氣動力對結構有限元厚度等設計變量的導數(shù)時，對結構剛度矩陣的影響進行了充分考慮。由于氣動彈性系統(tǒng)差分計算量極為龐大，因此本文給出了阻力系數(shù)對典型的設計變量(結構有限元設計變量主要對蒙皮與梁厚度)梯度與差分的對比，如圖15所示，橫坐標代表設計變量的編號，縱坐標代表目標函數(shù)對設計變量的梯度，耦合伴隨系統(tǒng)計算的梯度與差分(FD)幅值、趨勢基本一致，滿足多學科綜合設計要求,δ為擾動量。本文完成了氣動力對外形設計變量以及結構有限元設計變量兩類導數(shù)的驗證。結構應力(最常關注的是馮·米塞斯應力)對外形設計變量、結構設計變量的兩類導數(shù)，將在實現(xiàn)求解系統(tǒng)與有限元分析程序模塊緊耦合后進行校核，此時耦合伴隨系統(tǒng)右端項將更換為結構應力對多物理場狀態(tài)變量的變分形式，這是本文進一步的研究工作。

圖14 結構伴隨系統(tǒng)殘差收斂歷程Fig.14 Process of residual convergence of structural adjoint system

圖15 阻力系數(shù)對外形/結構設計變量的梯度對比Fig.15 Gradient of drag coefficients vs configuration/structural design variables

6 結 論

基于自主研發(fā)的大規(guī)模并行化結構網格RANS求解器PMB3D以及流固耦合模塊，開展了多學科耦合伴隨方程構造、求解方法的研究。

1) 基于流固耦合與流場求解共用模塊的不完全收斂解進行靜氣動彈性求解，具有較高的計算效率，相對于單次流場計算，計算量增加很少。

2) 耦合伴隨方程的第1項非對角線項可以通過流固耦合方式推導出其解析表達式。所采用的變形網格、參數(shù)化以及流固耦合方法在耦合伴隨系統(tǒng)推導中，決定了耦合系統(tǒng)本身的運算效率以及各項的具體表達形式。基于所采用的具體算法，本文給出了具體的推導結果，結合對原始變量變分以及鏈式求導，推導了虛功原理載荷傳遞方式下，流場伴隨方程交叉導數(shù)項的解析表達式；對變形網格混合算法求導問題采用有限差分近似簡化處理方式，降低了推導難度。

3) 結構伴隨方程的強迫項與梯度計算公式在表達形式上完全一致，僅在自變量擾動上不同，前者是結構位移，后者是設計變量，因此，文中結構伴隨方程強迫項計算與導數(shù)計算采用同一模塊來完成。文中進一步采用LDLT并行分解方法代替?zhèn)鹘y(tǒng)迭代方法，實現(xiàn)結構伴隨方程的高效方便求解，結構伴隨變量解趨于收斂僅需較少的耦合次數(shù)。

4) 延遲耦合伴隨系統(tǒng)能夠大大降低求解體系對內存的需求，同時也有利于程序的簡捷化、模塊化。

5) 文中采用的求解方法以及簡化處理方式所構建的離散伴隨系統(tǒng)效率較高，收斂穩(wěn)定性好，耦合系統(tǒng)的梯度求解精度滿足綜合設計需求。

本文的研究工作完成了氣動力對外形設計變量、結構有限元設計變量兩類導數(shù)的驗證，為進一步的多學科敏感性分析提供了研究基礎。下一步的研究內容為在實現(xiàn)求解系統(tǒng)與有限元分析程序模塊緊耦合后，求解結構應力對外形設計變量、結構設計變量的兩類導數(shù)。

致 謝

在本文的研究工作中，大連理工大學陳飆松教授、張盛副教授等在結構有限元建模、有限元分析等方面給予了有益的建議，空氣動力學國家重點實驗室孫巖助理研究員提供了標模計算網格，在此表示感謝！

參 考 文 獻

[1] MARTINS J R R A. A coupled-adjoint method for high-fidelity. Aero-structural optimization[M]. Stanford: Stanford University Press, 2002.

[2] MADER C A, KENWAY G K W, MARTINS J R R A. Towards high-fidelity aerostructural optimization using a coupled adjoint approach[C]∥12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference. Reston, VA: AIAA, 2008.

[3] KENWAY G K W, MARTINS J R R A. Multipoint high-fidelity aerostructural optimization of a transport aircraft configuration[J]. Journal of Aircraft, 2014, 51(1): 144-160.

[4] LEOVIRIYAKIT K, JAMESON A. Case studies in aero-structural wing planform and section optimization: AIAA-2004-5372[R]. Reston, VA: AIAA, 2004.

[5] ABU-ZURAYK M, BREZILLON J. Shape optimization using the aero-structural coupled adjoint approach for viscous flows[C]∥Evolutionary and Deterministic Methods For Design, Optimization and Control, 2011.

[6] MARCELET M, PETER J, CARRIER G. Sensitivity analysis of a strongly coupled aero-structural system using direct and adjoint methods: AIAA-2008-5863[R]. Reston, VA: AIAA, 2008.

[7] GHAZLANE I, CARRIER G, DUMONT A. Aerostructural adjoint method for flexible wing optimization: AIAA-2012-1924[R]. Reston, VA: AIAA, 2012.

[8] 左英桃， 高正紅， 詹浩.基于N-S方程和離散共軛方法的氣動設計方法研究[J].空氣動力學學報, 2009, 27(1):67-72.

ZUO Y T, GAO Z H, ZHAN H. Aerodynamic design method based on N-S equations and discrete adjoint approach[J]. Acta Aerodynamica Sinica, 2009, 27(1): 67-72 (in Chinese).

[9] 熊俊濤， 喬志德， 楊旭東， 等. 基于黏性伴隨方法的跨聲速機翼氣動優(yōu)化設計[J]. 航空學報, 2007, 28(2): 281-285.

XIONG J T, QIAO Z D, YANG X D, 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).

[10] 屈崑， 李記超， 蔡晉生. CFD數(shù)學模型的線性化方法及其應用[J]. 航空學報， 2015, 36(10): 3218-3227.

QU K, LI J C, CAI J S. Method of linearizing computational fluid dynamics model and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3218-3227 (in Chinese).

[11] 張朝磊， 厲海濤， 豐鎮(zhèn)平. 基于離散伴隨方法的透平葉柵氣動優(yōu)化[J]. 工程熱物理學報, 2012, 33(4): 47-50.

ZHANG C L, LI H T, FENG Z P. Aerodynamic opti-mization design of turbomachinery cascade based on discrete adjoint method[J]. Journal of Engineering Thermophysics, 2012, 33(4): 47-50 (in Chinese).

[12] 高宜勝， 伍貽兆， 夏健. 基于非結構網格離散型伴隨方法的翼型優(yōu)化[J]. 空氣動力學學報, 2013, 31(2): 244-249.

GAO Y S, WU Y Z, XIA J. A discrete adjoint-based approach for airfoil optimization on unstructured meshes[J]. Acta Aerodynamica Sinica, 2013, 31(2): 244-249 (in Chinese).

[13] 李彬， 鄧有奇， 唐靜， 等. 基于三維非結構混合網格的離散伴隨優(yōu)化方法[J]. 航空學報， 2014, 35(3): 674-686.

LI B, DENG Y Q, TANG J, et al. Discrete adjoint optimization method for 3D unstructured grid[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 674-686 (in Chinese).

[14] 張科施， 韓忠華， 李為吉， 等. 一種考慮氣動彈性的運輸機機翼多學科優(yōu)化方法[J]. 空氣動力學學報, 2008, 26(1): 1-7.

ZHANG K S, HAN Z H, LI W J, et al. A method of coupled aerodynamic/structural integration optimization for transport-wing design[J]. Acta Aerodynamica Sinica, 2008, 26(1): 1-7 (in Chinese).

[15] 劉剛, 肖中云, 王建濤, 等. 考慮約束的機載導彈導軌發(fā)射數(shù)值模擬[J]. 空氣動力學學報， 2015, 33(2): 192-197.

LIU G, XIAO Z Y, WANG J T, et al. Numerical simulation of missile air-launching process under rail slideway constraints[J]. Acta Aerodynamica Sinica, 2015, 33(2): 192-197 (in Chinese).

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

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 (in Chinese).

[17] ALLEN C B, RENDALL T C S. Unified approach to CFD-CSD interpolation and mesh motion using radial basis functions: AIAA-2007-3804[R]. Reston, VA: AIAA, 2007.

[18] 覃華， 徐燕子. 用LDLT并行分解優(yōu)化大規(guī)模SVM的訓練效率[J]. 計算機工程與應用, 2011, 47(12): 200-212.

QIN H， XU Y Z. Improving SVM’s learning efficiency by using matrix LDLT parallel decomposition[J]. Computer Engineering and Applications, 2011, 47(12): 200-212 (in Chinese).

[19] SPEKREIJSE S P, BOERSTOEL J W. An algorithm to check the topological validity of multi-block domain decompositions[C]∥Proceedings 6th International Conference on Numerical Grid Generation in Computational Field Simulations, 1998.

[20] BALLMANN J, BOUCKE A, CHEN B, et al. Aero-structural wind tunnel experiments with elastic wing models at high reynolds numbers (HIRENASD-ASDMAD): AIAA-2011-0882[R]. Reston, VA: AIAA, 2011.