飛行器氣動(dòng)/結(jié)構(gòu)多學(xué)科延遲耦合伴隨系統(tǒng)數(shù)值研究

黃江濤，周鑄,*，劉剛，高正紅，黃勇，王運(yùn)濤

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000 2.西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室， 西安 710072

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

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

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

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

1.1 并行化CFD代碼PMB3D

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

(1)

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

1.2 流固耦合代碼FSC3D

FSC3D代碼采用虛功原理與徑向基函數(shù)(RBF)插值技術(shù)[17]實(shí)現(xiàn)CFD邊界面元?dú)鈩?dòng)力向結(jié)構(gòu)點(diǎn)、結(jié)構(gòu)點(diǎn)位移向CFD物面格點(diǎn)高精度插值，基于RBF方法的插值技術(shù)數(shù)學(xué)模型表示為

(2)

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

本文對(duì)p(x)的選取采用多項(xiàng)式方法，即

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

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

采用LDLT方法[18]求解結(jié)構(gòu)靜力學(xué)方程：

Kd=F

(3)

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

(4)

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

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

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

2 流場(chǎng)伴隨方程簡(jiǎn)介

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

minI(w,X,D)

(5)

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

L=I+λTR

(6)

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

(7)

(8)

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

(9)

(10)

本文采用二階精度的中心格式、人工黏性以及SST兩方程湍流模型，手工推導(dǎo)構(gòu)造Navier-Stokes方程離散伴隨方程，并進(jìn)一步引入虛擬時(shí)間項(xiàng)：

(11)

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

Δλ*=0, Δλ=0

(12)

式中:λ*為隱式推進(jìn)的中間項(xiàng)。

黏性項(xiàng)采用薄層近似，曲線坐標(biāo)系下的黏性通量可以表達(dá)為

(13)

式中:各項(xiàng)表達(dá)式具體含義可參考文獻(xiàn)[16]。

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

3 考慮結(jié)構(gòu)靜力學(xué)方程的耦合伴隨方程

3.1 耦合伴隨方程

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

minI(w,X,D)

(14)

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

(15)

(16)

令式(16)第1項(xiàng)為0，獲取氣動(dòng)/結(jié)構(gòu)多學(xué)科耦合伴隨方程為

(17)

即

(18)

3.2 延遲耦合伴隨方程

(19)

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

其中，流場(chǎng)伴隨方程強(qiáng)迫項(xiàng)為

(20)

(21)

式中:wρ、wρ u、wρ v、wρ E和wρ w對(duì)應(yīng)守恒變量的展開(kāi)表達(dá)式。

由結(jié)構(gòu)靜力學(xué)方程(結(jié)構(gòu)節(jié)點(diǎn)個(gè)數(shù)為Ns),有

(22)

上述矩陣維度取決于結(jié)構(gòu)節(jié)點(diǎn)的維度以及CFD格心單元的維度，ψs的維度取決于有限元模型的網(wǎng)格類(lèi)型，對(duì)于板殼單元，ψs為6Ns×1維度的矢量，?Rs/?w的轉(zhuǎn)置矩陣As為5Na×6Ns矩陣。因此,結(jié)構(gòu)殘差對(duì)流場(chǎng)變量的導(dǎo)數(shù)?Rs/?w通過(guò)有限差分求出，解析表達(dá)式推導(dǎo)方式依賴于所采用的流固耦合方法。基于本文的虛功原理流固耦合方式，下面給出具體的推導(dǎo)。

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

(23)

(24)

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

(25)

代入耦合伴隨方程，可以獲得流場(chǎng)伴隨方程強(qiáng)迫項(xiàng)最終表達(dá)形式。

同理，結(jié)構(gòu)伴隨方程強(qiáng)迫項(xiàng)為

(26)

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

(27)

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

延遲伴隨方程第1項(xiàng)強(qiáng)迫項(xiàng)為大型矩陣與矢量相乘，將其具體表達(dá)式展開(kāi)，下標(biāo)表達(dá)式為

結(jié)構(gòu)伴隨方程強(qiáng)迫項(xiàng)的下標(biāo)表達(dá)式為

(28)

綜合強(qiáng)迫項(xiàng)的整理，可以推導(dǎo)出耦合延遲伴隨方程偽時(shí)間-殘差表達(dá)形式為

(29)

式(29)流場(chǎng)伴隨方程與結(jié)構(gòu)伴隨方程一步可以利用式(11)中的LU-SGS隱式迭代或雅克比迭代方法進(jìn)行求解。

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

分解出Lij和dii元素后，結(jié)構(gòu)伴隨方程求解可以通過(guò)三步掃描快速完成，即

(30)

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

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

對(duì)于延遲方程組式(19)中結(jié)構(gòu)伴隨方程的右端項(xiàng)，其含義是計(jì)算目標(biāo)函數(shù)對(duì)結(jié)構(gòu)位移的直接導(dǎo)數(shù)以及流場(chǎng)殘差對(duì)結(jié)構(gòu)位移導(dǎo)數(shù)與流場(chǎng)延遲伴隨變量的乘積，觀察單一流場(chǎng)學(xué)科目標(biāo)函數(shù)導(dǎo)數(shù)的計(jì)算公式：

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

(31)

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

圖3 流場(chǎng)單學(xué)科目標(biāo)函數(shù)導(dǎo)數(shù)計(jì)算模塊Fig.3 Objective function derivative module of flow field single discipline

圖4 結(jié)構(gòu)伴隨方程右端項(xiàng)計(jì)算模塊Fig.4 Computational module for right end item of structural adjoint equation

進(jìn)一步利用鏈?zhǔn)角髮?dǎo)展開(kāi)，即

(32)

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

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

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

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

(33)

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

對(duì)各項(xiàng)進(jìn)行進(jìn)一步展開(kāi)，有

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

4 學(xué)科分析耦合系統(tǒng)的驗(yàn)證

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

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

圖5 CFD網(wǎng)格分布與彈性變形Fig.5 CFD grid distribution and aeroelastic deformation

圖6 壓力分布與風(fēng)洞試驗(yàn)對(duì)比Fig.6 Pressure distribution compared with wind tunnel test

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

基于CRM (Common Research Model)標(biāo)準(zhǔn)算例構(gòu)型，進(jìn)行本文的方法梯度計(jì)算精度以及收斂特性驗(yàn)證。

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

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

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

圖7 CFD表面網(wǎng)格分布Fig.7 CFD gird distribution on wall

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

圖9 蒙皮結(jié)構(gòu)有限元網(wǎng)格Fig.9 Finite element grid of skin structure

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

圖11 氣動(dòng)彈性收斂歷程Fig.11 Process of aeroelastic convergence

圖12 典型耦合步長(zhǎng)彈性變形及局部視圖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

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

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

圖15 阻力系數(shù)對(duì)外形/結(jié)構(gòu)設(shè)計(jì)變量的梯度對(duì)比Fig.15 Gradient of drag coefficients vs configuration/structural design variables

6 結(jié) 論

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

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

2) 耦合伴隨方程的第1項(xiàng)非對(duì)角線項(xiàng)可以通過(guò)流固耦合方式推導(dǎo)出其解析表達(dá)式。所采用的變形網(wǎng)格、參數(shù)化以及流固耦合方法在耦合伴隨系統(tǒng)推導(dǎo)中，決定了耦合系統(tǒng)本身的運(yùn)算效率以及各項(xiàng)的具體表達(dá)形式?；谒捎玫木唧w算法，本文給出了具體的推導(dǎo)結(jié)果，結(jié)合對(duì)原始變量變分以及鏈?zhǔn)角髮?dǎo)，推導(dǎo)了虛功原理載荷傳遞方式下，流場(chǎng)伴隨方程交叉導(dǎo)數(shù)項(xiàng)的解析表達(dá)式；對(duì)變形網(wǎng)格混合算法求導(dǎo)問(wèn)題采用有限差分近似簡(jiǎn)化處理方式，降低了推導(dǎo)難度。

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

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

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

本文的研究工作完成了氣動(dòng)力對(duì)外形設(shè)計(jì)變量、結(jié)構(gòu)有限元設(shè)計(jì)變量?jī)深?lèi)導(dǎo)數(shù)的驗(yàn)證，為進(jìn)一步的多學(xué)科敏感性分析提供了研究基礎(chǔ)。下一步的研究?jī)?nèi)容為在實(shí)現(xiàn)求解系統(tǒng)與有限元分析程序模塊緊耦合后，求解結(jié)構(gòu)應(yīng)力對(duì)外形設(shè)計(jì)變量、結(jié)構(gòu)設(shè)計(jì)變量的兩類(lèi)導(dǎo)數(shù)。

致 謝

在本文的研究工作中，大連理工大學(xué)陳飆松教授、張盛副教授等在結(jié)構(gòu)有限元建模、有限元分析等方面給予了有益的建議，空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室孫巖助理研究員提供了標(biāo)模計(jì)算網(wǎng)格，在此表示感謝！

參 考 文 獻(xiàn)

[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方程和離散共軛方法的氣動(dòng)設(shè)計(jì)方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào), 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] 熊俊濤， 喬志德， 楊旭東， 等. 基于黏性伴隨方法的跨聲速機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 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ù)學(xué)模型的線性化方法及其應(yīng)用[J]. 航空學(xué)報(bào)， 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)平. 基于離散伴隨方法的透平葉柵氣動(dòng)優(yōu)化[J]. 工程熱物理學(xué)報(bào), 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] 高宜勝， 伍貽兆， 夏健. 基于非結(jié)構(gòu)網(wǎng)格離散型伴隨方法的翼型優(yōu)化[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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] 李彬， 鄧有奇， 唐靜， 等. 基于三維非結(jié)構(gòu)混合網(wǎng)格的離散伴隨優(yōu)化方法[J]. 航空學(xué)報(bào)， 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] 張科施， 韓忠華， 李為吉， 等. 一種考慮氣動(dòng)彈性的運(yùn)輸機(jī)機(jī)翼多學(xué)科優(yōu)化方法[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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] 劉剛, 肖中云, 王建濤, 等. 考慮約束的機(jī)載導(dǎo)彈導(dǎo)軌發(fā)射數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)， 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] 黃江濤， 劉剛， 周鑄, 等. 基于離散伴隨方程求解梯度信息的若干問(wèn)題研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào)， 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的訓(xùn)練效率[J]. 計(jì)算機(jī)工程與應(yīng)用, 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.