王晨, 白俊強(qiáng), Jan S Hesthaven, 邱亞松, 喬磊, 韓嘯

(1.西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2.洛桑聯(lián)邦理工(EPFL), 瑞士 洛桑 CH-1015)

在航空航天領(lǐng)域，對(duì)于復(fù)雜的帶參數(shù)非定常非線性問(wèn)題，基于本征正交分解方法(POD)的降階模型(ROM)[1-2]被廣泛采用以提高計(jì)算效率，其思路是通過(guò)POD對(duì)復(fù)雜非線性解空間構(gòu)造一組最優(yōu)降階基以保留流場(chǎng)的主要特征，再將待求解空間投影到由這組基張成的線性子空間中以快速預(yù)測(cè)未知流場(chǎng)的低階模型。常用的投影方法包括Galerkin投影等，但對(duì)于處理復(fù)雜的帶參數(shù)非線性問(wèn)題時(shí)，采用投影方法給投影系數(shù)求解帶來(lái)不便，為此，目前有一批研究學(xué)者考慮采用KrigingRBF等代理模型來(lái)預(yù)測(cè)降階模型中各POD基投影系數(shù)的降階模型形式，如張偉偉等[3]則采用POD與徑向基函數(shù)法(radial basis function,RBF)的結(jié)合預(yù)測(cè)翼型的氣動(dòng)響應(yīng)和極限環(huán)振蕩，邱亞松等[4]采用Kriging代理模型預(yù)測(cè)POD基投影系數(shù)，對(duì)由10個(gè)參數(shù)控制的翼型擾流流場(chǎng)進(jìn)行了預(yù)測(cè)。本文考慮采用POD結(jié)合Kriging代理模型的形式構(gòu)造降階模型，預(yù)測(cè)帶參數(shù)的典型流動(dòng)問(wèn)題——串聯(lián)雙圓柱構(gòu)型。

在數(shù)值方法和流動(dòng)機(jī)理的研究中，圓柱擾流是個(gè)經(jīng)典的非定常算例。對(duì)于計(jì)算狀態(tài)確定的等直徑串聯(lián)雙圓柱，距徑比L/D(圓柱柱心間距與直徑的比值，L為圓柱柱心間距；D為圓柱直徑)是唯一影響流動(dòng)形態(tài)的參數(shù)。美國(guó)Langley實(shí)驗(yàn)室(LaRC)提供了L/D=3.7、L/D=1.435下basic aerodynamic research tunnel(BART)[5]和quiet flow facility(QFF)[6]2組風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)。Zdravkovich等[7-8]通過(guò)研究，總結(jié)了不同距徑比下雙圓柱流場(chǎng)的流動(dòng)形態(tài)規(guī)律。由于雙圓柱流動(dòng)具備強(qiáng)烈的非定常效應(yīng)，又有對(duì)流動(dòng)形態(tài)起著決定性作用的幾何參數(shù)，并有實(shí)驗(yàn)數(shù)據(jù)作支撐，是降階模型理想的驗(yàn)證算例。本文試圖通過(guò)降階模型對(duì)變距徑比串聯(lián)雙圓柱流場(chǎng)進(jìn)行預(yù)測(cè)。

但在采用POD提取降階基之前，降階模型要求提供樣本流場(chǎng)解所組成的快照數(shù)據(jù)，為此需確定雙圓柱流動(dòng)的數(shù)值模擬方法計(jì)算樣本數(shù)據(jù)。對(duì)于串聯(lián)雙圓柱構(gòu)型，很多學(xué)者對(duì)BART和QFF的實(shí)驗(yàn)工況采用不同的數(shù)值模擬方法進(jìn)行了相關(guān)研究：Khorrami等[9-10]采用基于SST湍流模型的非定常雷諾平均方法(URANS)；Wessels等[11]采用格子波爾茲曼法分別對(duì)LaRC風(fēng)洞實(shí)驗(yàn)工況進(jìn)行非定常數(shù)值仿真。而相對(duì)于URANS方法，能解析非定常脈動(dòng)的LES、混合RANSLES方法在模擬雙圓柱流動(dòng)時(shí)將更為精確，如白俊強(qiáng)等[12]通過(guò)采用尺度自適應(yīng)模型(SAS)[13]對(duì)L/D=3.7的臨界工況進(jìn)行模擬，說(shuō)明了SAS不僅能捕捉瞬時(shí)脈動(dòng)，在時(shí)均結(jié)果的模擬精度上也優(yōu)于URANS方法。本文選擇SAS方法提供快照數(shù)據(jù)。

針對(duì)串聯(lián)雙圓柱流動(dòng)，本文采用POD提取主要的降階基模態(tài)以構(gòu)造線性降階子空間，再利用Kriging代理模型插值投影系數(shù)來(lái)封閉降階模型以預(yù)測(cè)未知流場(chǎng)。通過(guò)與風(fēng)洞實(shí)驗(yàn)及數(shù)值求解計(jì)算結(jié)果作對(duì)比，評(píng)估降階模型在變參大分離類問(wèn)題上的預(yù)測(cè)精度與效率。為保證快照的數(shù)值精度，快照數(shù)據(jù)由非定常數(shù)值求解器采集，模型選擇的是基于SST[14]的尺度自適應(yīng)模型(SAS)。

1 降階基模型建立

對(duì)于帶參數(shù)定常問(wèn)題uθ(x)(x∈Ω),參數(shù)向量θ,系統(tǒng)提供N個(gè)樣本uθ(i)(x)(i=1,2,…N),降階基模型通過(guò)從樣本中提取一組正交降階基Φk(k=1,2,…,K)來(lái)預(yù)測(cè)結(jié)果:

θ(x)=u0(x)+∑Kk=1αk(θ)·Φk(x)

(1)

(Φi,Φj)=0i≠j

1i=j(i,j=1,2,…K)

(2)

αk(θ)是待求的投影系數(shù)。u0(x)通常取為快照的平均u0(x)=1N∑Ni=1uθ(i)(x),這一項(xiàng)的添加能夠顯著提高降階基模型的預(yù)測(cè)精度。

同樣對(duì)于非定常非線性問(wèn)題u(x,t),降階基模型給出類似的預(yù)測(cè)公式:

(x,t)=u0(x)+∑Kk=1αk(t)·Φk(x)

(3)

1.1 降階基提取

(1)、(3)式中的降階基Φk(x)可由不同的方法確定,本文選用本征正交分解(POD)方法提取Φk(x)。

POD方法提取的降階基Φk(x)需滿足快照在降階子空間上的投影達(dá)到最大,即:

A·v=λv(5)

矩陣A是對(duì)稱非負(fù)定矩陣,具有N個(gè)特征值λk和相應(yīng)的特征向量vk,與之對(duì)應(yīng)的基模態(tài)Φk(x)可通過(guò)k(x)=∑Ni=1vk(i)·uθ(i)(x)再對(duì)k(x)單位化后求得。為節(jié)省計(jì)算量,需要從這N個(gè)模態(tài)中選擇最主要的模態(tài)來(lái)張成降階子空間。基于j(x),j(x)〉dx=∑Nk=1λk(u′(·)=u(·)-u0(·)),Sirovich[15]提出了基于特征值λk的廣義能法來(lái)選擇主要模態(tài):將特征值降序排列,定義能量百分比I(k)=∑ki=1λi∑Ni=1λi,給定能量閾值ε∈(0,1),通常取0.99,則K可通過(guò)式K=argmin(I(k)I(k)≥ε)確定,前K個(gè)特征值對(duì)應(yīng)的模態(tài)就是最終獲得的POD模態(tài)Φk(x)(k=1,2,…,K)。

1.2 基系數(shù)預(yù)測(cè)

Kriging代理模型[16]是通過(guò)對(duì)樣本的訓(xùn)練建立輸入、輸出函數(shù)關(guān)系,對(duì)輸出變量給出最佳線性無(wú)偏估計(jì)的插值方法。在本文中,Kriging模型并非直接用來(lái)預(yù)測(cè)狀態(tài)變量(若將uθ)(x)(x∈Ω)直接作為輸出,就需對(duì)每個(gè)網(wǎng)格上的狀態(tài)變量建立Kriging模型,那么代理模型將會(huì)失去其原有的效率優(yōu)勢(shì)),而是將提取的降階基系數(shù)αk(θ)作為輸出,這樣通過(guò)Kriging模型提供投影系數(shù)αk(θ),然后再通過(guò)降階模型預(yù)測(cè)公式(1)、(3)預(yù)測(cè)待求流場(chǎng)。對(duì)單個(gè)狀態(tài)變量來(lái)說(shuō),POD提取了K個(gè)降階基,就需建立K個(gè)代理模型。在用POD提取基模態(tài)后,N個(gè)樣本的投影系數(shù)可以通過(guò)下式得到:

αk(θ(i))=Φk(x)·(uθ(i)(x)-u0(x))

i=1,2,…,N(6)

以此作為代理模型的樣本,每個(gè)Kriging代理模型的構(gòu)造和預(yù)測(cè)過(guò)程如下:

設(shè)輸出α,輸入為p維向量θ,由1.1節(jié)得到對(duì)應(yīng)的N組樣本,Kriging模型對(duì)輸入輸出建立了如下形式的函數(shù)關(guān)系:

α(θ)=∑pi=1βifi(θ)+Z(θ)

(7)

式中，f(θ)為回歸模型,Z(θ)為期望為0、方差為σ2、協(xié)方差為cov(Z(θ(i)),Z(θ(j)))=σ2R(θ(i),θ(j))的正態(tài)分布。θ(i)為第i個(gè)樣本,相關(guān)函數(shù)R表征著樣本兩兩之間在θ空間內(nèi)的相關(guān)程度,其計(jì)算公式如下:

R(θ(i),θ(j))=e(-d(θ(i),θ(j)))

min:ψ(wT)=|R|1/pσ2

(θ)=f(θ)·+(R*(θ))T·R-1·(A-F)

R*(θ)=[R(θ,θ(1)),R(θ,θ(2)),…,R(θ,θ(N))]T

(10)

通過(guò)K個(gè)代理模型的建立,就可以確定參數(shù)域內(nèi)任一點(diǎn)降階模型的系數(shù)αk(θ),再由(1)、(3)式即可預(yù)測(cè)出對(duì)應(yīng)的流場(chǎng)數(shù)據(jù)。

2 串聯(lián)雙圓柱流動(dòng)簡(jiǎn)介及數(shù)值模擬方法

2.1 串聯(lián)雙圓柱流動(dòng)簡(jiǎn)介

計(jì)算工況采用NASA雙圓柱風(fēng)洞實(shí)驗(yàn)工況(見(jiàn)圖1):來(lái)流馬赫數(shù)0.128 5(V=44 m/s),基于圓柱直徑的雷諾數(shù)為1.66×105,圓柱直徑D=0.057 1 m?？紤]到三維非定常計(jì)算過(guò)于耗時(shí),本文采用二維非定常計(jì)算。圓柱表面采用無(wú)滑移邊界條件,遠(yuǎn)場(chǎng)采用無(wú)反射邊界條件。物理時(shí)間步長(zhǎng)取Δt=1.0×10-5s。

圖1 串聯(lián)雙圓柱幾何示意圖

2.2 非定常數(shù)值模擬方法

本文采用雙時(shí)間步,空間離散格式為二階中心差分的有限體積求解器對(duì)雙圓柱流場(chǎng)進(jìn)行非定常模擬。

因?yàn)殡p圓柱的流場(chǎng)中存在明顯的渦旋結(jié)構(gòu),而非定常雷諾平均方法(URANS)對(duì)分離流動(dòng)的預(yù)測(cè)精度低,本文采用Menter在k-ωSST湍流模型的基礎(chǔ)上發(fā)展的一種混合RANSLES模型——尺度自適應(yīng)模型(SAS)方法來(lái)保證湍流預(yù)測(cè)的精度。流動(dòng)方程是濾波后的Navier-Stokes方程:

?ρ?t+?(ρU)?x+?(ρV)?y=0

?(ρU)?t+?(ρU2)?x+?(ρUV)?y=-?p?x+

(μ+μt)?2U?x2+?2U?y2

?(ρV)?t+?(ρUV)?x+?(ρV2)?y=-?p?y+

(μ+μt)?2V?x2+?2V?y2

?ρe+U2+V22?t+·ρe+U2+V22U=-

(11)

式中，亞格子黏性系數(shù)μt,表征小尺度脈動(dòng)對(duì)于大尺度脈動(dòng)的作用,其求解仍采用了k-ωSST求解渦粘性系數(shù)的表達(dá)形式:

μt=minρkω,a1ρkΩF2

(12)

?(ρk)?t+ui?(ρk)?xi=Pk-

βkρkω+??xiμ+μtσk?k?xi

?(ρω)?t+ui?(ρω)?xi=CωPω+QSAS-

βωρω2+??xiμ+μtσω?ω?xi+

2ρ(1-F1)1σω21ω?k?xi?ω?xi

(13)

但在求解湍動(dòng)能k和湍流比耗散率ω時(shí),混合RANSLES通過(guò)引入比湍流模型中積分尺度更小的特征尺度,以保證模型具備識(shí)別非定常脈動(dòng)的能力。其中DES類方法利用網(wǎng)格尺度增大湍動(dòng)能(k)方程的耗散,而本文選用的SAS方法則是通過(guò)在比耗散率(ω)方程中添加由馮卡門尺度確定的源項(xiàng)QSAS來(lái)解析流動(dòng)的脈動(dòng)效應(yīng),保證模型具備捕捉非定?，F(xiàn)象的能力。

SAS添加源項(xiàng):

QSAS=ρFSASmaxζ2S2LLvk-2kσφFSST-SAS,0(14)

馮卡門尺度Lvk=κU′U″,表征流場(chǎng)當(dāng)?shù)厮俣鹊囊浑A導(dǎo)和二階導(dǎo)的比值,它與網(wǎng)格尺度無(wú)關(guān)。用馮卡門尺度取代網(wǎng)格尺度來(lái)減小模型識(shí)別的流動(dòng)特征尺度,使得SAS對(duì)網(wǎng)格的依賴程度要低于DES類方法。

2.3 變參數(shù)計(jì)算網(wǎng)格生成

對(duì)于不同距徑比的雙圓柱構(gòu)型,在進(jìn)行數(shù)值計(jì)算時(shí),本文均采用相同的計(jì)算網(wǎng)格拓?fù)?。網(wǎng)格采用多塊結(jié)構(gòu)網(wǎng)格,總計(jì)單元數(shù)11.2萬(wàn),如圖2所示:邊界層內(nèi)部仍采用RANS的網(wǎng)格要求,即物面法向第一層網(wǎng)格滿足保證y+<1,附面層內(nèi)網(wǎng)格物面法向增長(zhǎng)率為1.1;附面層網(wǎng)格外部進(jìn)行適當(dāng)加密以滿足非定常模擬精度。

圖2 計(jì)算網(wǎng)格示意圖

本文先以圖2的網(wǎng)格拓?fù)鋵?duì)L/D=3.7的工況建立計(jì)算網(wǎng)格,當(dāng)距徑比變化時(shí),通過(guò)程序生成新的雙圓柱幾何,再用無(wú)限代數(shù)插值(TFI)[17]動(dòng)網(wǎng)格技術(shù)生成新構(gòu)型的計(jì)算網(wǎng)格。

在進(jìn)行計(jì)算時(shí),CFD求解器首先將物理坐標(biāo)系(oxy)轉(zhuǎn)換到單位正交計(jì)算坐標(biāo)系(o′x′y′)中進(jìn)行計(jì)算,如圖3所示。這樣對(duì)于不同距徑比構(gòu)型采用相同網(wǎng)格拓?fù)浜?盡管在物理坐標(biāo)系下各網(wǎng)格點(diǎn)的空間坐標(biāo)不同,但其在計(jì)算坐標(biāo)系下坐標(biāo)是完全相同的,本文對(duì)變距徑比構(gòu)型采用(1)式、(3)式建立降階模型時(shí)采用的正是計(jì)算坐標(biāo)系。采用相同的網(wǎng)格拓?fù)?和計(jì)算坐標(biāo)系的使用使得降階模型可以處理變幾何參數(shù)類問(wèn)題。

圖3 坐標(biāo)系轉(zhuǎn)換示意圖

3 降階模型在串聯(lián)雙圓柱流動(dòng)中的應(yīng)用與預(yù)測(cè)

3.1 變參數(shù)時(shí)均場(chǎng)模擬

作為串聯(lián)雙圓柱構(gòu)型重要的流場(chǎng)參數(shù),距徑比L/D的數(shù)值直接決定著雙圓柱流場(chǎng)形態(tài),為評(píng)估降階模型模擬變參數(shù)流場(chǎng)的能力,本文試圖通過(guò)一系列L/D的時(shí)均流場(chǎng)構(gòu)造降階模型,預(yù)測(cè)L/D=3.7.1.435這2個(gè)有實(shí)驗(yàn)數(shù)據(jù)的臨界距比的流場(chǎng),以評(píng)估模型對(duì)于預(yù)測(cè)變參類問(wèn)題的精度與效率。

本文采用2.2節(jié)的非定常數(shù)值模擬方法和2.3節(jié)的網(wǎng)格拓?fù)涮峁?0組不同L/D的時(shí)均流場(chǎng)作為快照數(shù)據(jù):{1.2,1.25,1.3,1.35,1.55,1.7,1.85,2,2.25,2.5,2.75,3,3.25,3.5,3.6,3.85,4,4.25,4.5,4.75,5}。采用POD對(duì)4個(gè)狀態(tài)變量(ρ,u,v,p)的快照分別提取降階基,再通過(guò)廣義能法選取主要模態(tài):廣義能百分比隨基數(shù)目的變化曲線如圖4所示,圖中對(duì)于4個(gè)狀態(tài)變量,前5個(gè)基模態(tài)占據(jù)了幾乎所有的廣義能;給定99.5%的能量閾值,本文最終為ρ,u,v,p分別選取5,5,4,4個(gè)降階基。然后對(duì)這5+5+4+4=18個(gè)POD基的權(quán)重系數(shù)分別建立Kriging代理模型來(lái)預(yù)測(cè)POD基系數(shù),這樣對(duì)于任意L/D∈[1.25,5],結(jié)合降階基和代理模型,就能預(yù)測(cè)出L/D在1.25到5范圍內(nèi)任一距徑比的時(shí)均流場(chǎng)。

圖4 能量百分比隨基模態(tài)數(shù)目的變化曲線

1) 預(yù)測(cè)精度評(píng)估

本文重點(diǎn)考慮L/D=3.7和1.435的預(yù)測(cè)精度。

L/D=3.7和1.435時(shí)風(fēng)洞實(shí)驗(yàn)、數(shù)值求解器、降階模型三者得到的空間流線分布如圖5所示。

圖5 時(shí)均流線對(duì)比圖

其中實(shí)驗(yàn)流線圖來(lái)源于文獻(xiàn)[10]。L/D=3.7時(shí),數(shù)值模擬和降階基預(yù)測(cè)的下游圓柱空間背風(fēng)面的分離區(qū)域大于實(shí)驗(yàn)所探測(cè)的;L/D=1.435時(shí),實(shí)驗(yàn)數(shù)據(jù)中兩圓柱間是個(gè)整體的分離渦形式,而數(shù)值結(jié)果則在圓柱間預(yù)測(cè)出上下2個(gè)相互接觸的漩渦結(jié)構(gòu)。盡管與實(shí)驗(yàn)的時(shí)均流線存在差異,降階模型與求解器的流線包括圓柱背風(fēng)面分離區(qū)的大小和形態(tài)幾乎完全一致,這說(shuō)明時(shí)均流場(chǎng)的差異主要是由于數(shù)值求解方法所致,而并非降階模型預(yù)測(cè)誤差所造成。

為進(jìn)一步檢驗(yàn)降階模型的預(yù)測(cè)精度,考察降階模型預(yù)測(cè)的時(shí)均流場(chǎng)和數(shù)值求解之間的誤差分布。定義降階模型預(yù)測(cè)值與數(shù)值求解器計(jì)算值的絕對(duì)誤差Δ(x)=‖θ(x)-uθ(x)‖,相對(duì)誤差Δr(x)=‖θ(x)-uθ(x)‖uθ(x)。圖6給出ρ,p在L/D=3.7,1.453時(shí)的相對(duì)誤差分布,考慮到距徑比對(duì)遠(yuǎn)離物面處流動(dòng)的影響小,此時(shí)降階模型誤差小,故這里僅提供雙圓柱附近流場(chǎng)的誤差分布。在整個(gè)考察區(qū)域內(nèi)降階模型預(yù)測(cè)的相對(duì)誤差處在10-3量級(jí),說(shuō)明降階模型對(duì)ρ,p實(shí)現(xiàn)了精確的預(yù)測(cè)。至于速度場(chǎng),在本算例中,圓柱的背風(fēng)面是流動(dòng)滯止區(qū),在這個(gè)區(qū)域內(nèi)的速度相對(duì)其他區(qū)域數(shù)值很小,使得降階模型在該區(qū)域某些地方相對(duì)誤差過(guò)大,故本文考慮速度場(chǎng)的絕對(duì)誤差分布,如圖7所示。圖中流場(chǎng)各處速度的絕對(duì)誤差始終處小于10-2,而在圖5中降階模型又能較精確地模擬兩工況的時(shí)均流線,因此本文得出降階模型對(duì)速度場(chǎng)、密度和壓強(qiáng)均實(shí)現(xiàn)了全流場(chǎng)精確預(yù)測(cè)的結(jié)論。

物體表面壓強(qiáng)分布是氣動(dòng)計(jì)算與預(yù)測(cè)的關(guān)注點(diǎn)之一,降階模型預(yù)測(cè)的圓柱表面壓力分布如圖8、圖9所示。它與求解器的計(jì)算結(jié)果除在駐點(diǎn)處外幾乎完全重合,在L/D=3.7預(yù)測(cè)的兩圓柱壓力分布及L/D=1.435時(shí)上游圓柱的壓強(qiáng)系數(shù)分布與實(shí)驗(yàn)也保持一致,但L/D=1.435下游圓柱與實(shí)驗(yàn)則相差較大,這和未能模擬出實(shí)驗(yàn)中此距徑比下兩圓柱間的分離形態(tài)相吻合。這種差距一方面是因?yàn)閿?shù)值模擬方法的精度仍有待提高,另一方面是由于數(shù)值計(jì)算時(shí)沒(méi)有考慮風(fēng)洞實(shí)驗(yàn)的某些條件和外在因素所致。

圖6 時(shí)均流場(chǎng) 圖7 時(shí)均流場(chǎng)速度絕對(duì)誤差分布圖

圖8 L/D=3.7圓柱表面壓力分布對(duì)比 圖9 L/D=1.435圓柱表面壓力分布對(duì)比

通過(guò)時(shí)均流線、誤差分布、表面壓力分布的對(duì)比,本文的降階模型在對(duì)20個(gè)樣本提取降階基后,能精確地模擬L/D=3.7和L/D=1.435的時(shí)均流場(chǎng)。

2) 計(jì)算效率評(píng)估

除預(yù)測(cè)精度之外,計(jì)算效率是評(píng)估降階模型性能優(yōu)劣的另一個(gè)重要因素,降階模型各階段計(jì)算耗時(shí)如表1所示:

①降階模型對(duì)20個(gè)樣本生成降階子空間耗時(shí)22.42 s,對(duì)POD基系數(shù)建立代理模型時(shí)花費(fèi)1 452 s,這兩部分構(gòu)成了整個(gè)訓(xùn)練階段;

②在降階模型建立后,對(duì)于單點(diǎn)流場(chǎng),降階模型僅需2.3 s就能較精確地預(yù)測(cè)出時(shí)均結(jié)果。而采用求解器計(jì)算時(shí)均結(jié)果至少需要兩萬(wàn)非定常時(shí)間步,需耗時(shí)19 382 s。因此,本文為雙圓柱流場(chǎng)構(gòu)造的降階模型利用訓(xùn)練階段的基函數(shù)和代理模型能快速又精確地預(yù)測(cè)變距徑比時(shí)均流場(chǎng),其在計(jì)算效率上帶來(lái)的收益非?？捎^。

表1 計(jì)算效率評(píng)估

3) 時(shí)均流場(chǎng)變化趨勢(shì)分析

在L/D∈[1.25,5]范圍內(nèi)均勻選取750個(gè)測(cè)試點(diǎn),用建立的降階模型提取時(shí)均流場(chǎng)后,本文對(duì)変距徑比串聯(lián)雙圓柱時(shí)均流場(chǎng)得出以下規(guī)律(這750個(gè)測(cè)試點(diǎn)的網(wǎng)格數(shù)據(jù)是通過(guò)對(duì)空間坐標(biāo)另外建立POD-Kriging降階模型得到的):

①在L/D∈[1.25,2.815]范圍內(nèi),上游圓柱背風(fēng)面的分離渦占據(jù)著圓柱間的整個(gè)區(qū)域,如圖10a)所示;

②在L/D∈(2.815,5]范圍內(nèi),上游圓柱背風(fēng)面的分離渦僅分布在兩圓柱間的部分區(qū)域,其與下游圓柱背風(fēng)面并未接觸,如圖10b)所示。

圖10 時(shí)均流場(chǎng)形態(tài)

3.2 定參非定常模擬

3.1節(jié)對(duì)變距徑比時(shí)均流場(chǎng)進(jìn)行了分析,在本小節(jié)中,針對(duì)L/D=3.7這個(gè)臨界工況的非定常流場(chǎng)建立降階基,以考察降階模型對(duì)非定常問(wèn)題的預(yù)測(cè)精度。對(duì)該工況計(jì)算5萬(wàn)非定常時(shí)間步后,令此時(shí)t=0 s,再在t∈[0,0.5]范圍內(nèi)均勻選取501個(gè)時(shí)刻的瞬時(shí)流場(chǎng)作為快照數(shù)據(jù)。采用POD對(duì)各狀態(tài)變量(ρ,u,v,p)的快照分別提取降階基,能量百分比隨POD基數(shù)的變化曲線如圖11所示,利用廣義能法為ρ,u,v,p分別選取(15,15,11,13)個(gè)降階基。

圖11 能量百分比隨基模態(tài)數(shù)目的變化曲線

取t=0.101 5 s、0.206 5 s、0.336 5 s、0.473 5 s這4個(gè)時(shí)刻作為測(cè)試點(diǎn),顯然這4個(gè)測(cè)試點(diǎn)均不屬于樣本集。圖12為速度的絕對(duì)誤差分布,在考察的4個(gè)時(shí)刻降階模型預(yù)測(cè)的絕對(duì)誤差始終保持在10-3量級(jí),說(shuō)明了本文建立的降階模型在t∈[0,0.5]范圍內(nèi)對(duì)瞬時(shí)速度場(chǎng)具有較高的預(yù)測(cè)精度。圖13為測(cè)試時(shí)刻的瞬時(shí)渦量圖,左圖中本文采用的數(shù)值模擬方法能在一定程度上反映出捕捉流場(chǎng)的非定常特征,而對(duì)于降階模型來(lái)說(shuō),非定常效應(yīng)往往難以精確模擬,右圖中本文采用的降階模型預(yù)測(cè)的渦量分布與左圖數(shù)值計(jì)算的幾乎完全一致,證明了其對(duì)非線性、非定常問(wèn)題具有較高的預(yù)測(cè)精度。

圖12 測(cè)試點(diǎn)速度絕對(duì)分布圖13 瞬時(shí)渦量圖

4 結(jié) 論

本文通過(guò)采用POD提取主要的降階基模態(tài)，對(duì)流場(chǎng)解構(gòu)造了線性降階子空間，再利用Kriging代理模型插值投影系數(shù)來(lái)封閉降階模型，對(duì)變距徑比串聯(lián)雙圓柱流場(chǎng)進(jìn)行預(yù)測(cè)，得到以下結(jié)論：

1)在與實(shí)驗(yàn)對(duì)比后，本文采取的非定常求解方法能較精確預(yù)測(cè)串聯(lián)雙圓柱非定常流動(dòng)；

2)對(duì)于變距徑比構(gòu)型，降階模型能較精確預(yù)測(cè)出非樣本點(diǎn)的時(shí)均流場(chǎng)，相對(duì)于全階模型計(jì)算效率提高69 470倍，其與實(shí)驗(yàn)數(shù)據(jù)的誤差主要源于數(shù)值求解的誤差；

3)通過(guò)建立的降階模型，本文總結(jié)了雙圓柱時(shí)均流場(chǎng)隨距徑比的變化趨勢(shì)；

4)在特定工況下，POD-Kriging降階模型同樣能精確預(yù)測(cè)給定時(shí)刻的瞬時(shí)流場(chǎng)。

考慮到傳統(tǒng)的POD-Galerkin方法在多參數(shù)問(wèn)題中構(gòu)造不便，本文并未給出其與傳統(tǒng)的POD-Galerkin方法的對(duì)比而只給出了其與全階模型精度及計(jì)算效率的對(duì)比。另外，由于三維非定常計(jì)算非常耗時(shí)，為保證取樣效率以印證降階模型對(duì)變參串聯(lián)雙圓柱流場(chǎng)的預(yù)測(cè)精度與效率，本文采用二維非定常計(jì)算，但由于忽略展向流動(dòng)，這在一定程度勢(shì)必會(huì)損失預(yù)測(cè)精度，后續(xù)工作將會(huì)針對(duì)三維數(shù)值模擬結(jié)果建立適應(yīng)于帶參數(shù)非定常問(wèn)題的降階基模型，此時(shí)相對(duì)于二維流場(chǎng)預(yù)測(cè)，降階模型在計(jì)算效率上的收益預(yù)計(jì)將會(huì)更為可觀。

針對(duì)變參問(wèn)題的降階模型在實(shí)際應(yīng)用如氣動(dòng)設(shè)計(jì)、流固耦合等方向中將擁有比較廣泛和誘人的應(yīng)用前景，后續(xù)研究將集中在高維參數(shù)空間的強(qiáng)非線性非定常問(wèn)題中提高降階模型魯棒性和自適應(yīng)取樣這兩項(xiàng)工作。

參考文獻(xiàn)：

[1] Thomas J P, Dowell E H, Hall K C. Three-Dimensional Transonic Aeroelasticity Using Proper Orthogonal Decomposition Based Reduced Order Models[J]. Journal of Aircraft, 2003, 40(3): 544-551

[2] Du J, Fang F, Pain C C, et al. POD Reduced-Order Unstructured Mesh Modeling Applied To 2D and 3D Fluid Flow[J]. Computers & Mathematics with Applications, 2013, 65(3): 362-379

[3] Zhang W, Kou J, Wang Z. Nonlinear Aerodynamic Reduced-Order Model for Limit-Cycle Oscillation and Flutter[J]. AIAA Journal, 2016, 54(10): 3304-3311

[4] Qiu Yasong, Bai Junqiang. Stationary Flow Fields Prediction of Variable Physical Domain Based on Proper Orthogonal Decomposition and Kriging Surrogate Model[J]. Chinese Journal of Aeronautics, 2015, 18(1): 44-56

[5] Jenkins L N, Khorrami M R, Choudhari M M, et al. Characterization of Unsteady Flow Structures around Tandem Cylinders for Component Interaction Studies in Airframe Noise[R]. AIAA-2005-2812

[6] Jenkins L N, Neuhart D H, McGinley C B, et al. Measurements of Unsteady Wake Interference between Tandem Cylinders[R]. AIAA-2006-3202

[7] Zdravkovich M M. Flow around Circular Cylinders Volume 2: Applications[M]. Oxford University Press, 1997

[8] Zdravkovich M M. Review of Flow Interference between Two Circular Cylinders in Various Arrangements[J]. Journal of Fluids Engineering, 1977, 99(4): 618

[9] Khorrami M R, Choudhari M M, Lockard D P, et al. Unsteady Flowfield around Tandem Cylinders as Prototype Component Interaction in Airframe Noise[J]. AIAA Journal, 2007, 45(8): 1930-1941

[10] Lockard D P, Choudhari M M, Khorrami M R, et al. Aeroacosutic Simulations of Tandem Cylinders with Subcritical Spacing[R]. AIAA-2008-2862

[11] Brès G, Wessels M, Noelting S. Tandem Cylinder Noise Predictions Using Lattice Boltzmann and Ffowcs Williams-Hawkings Methods[R]. AIAA-2010-3791

[12] 白俊強(qiáng),王晨,張揚(yáng). 一種基于馮卡門尺度的湍流模式在模擬穩(wěn)態(tài)和非穩(wěn)態(tài)流動(dòng)問(wèn)題中的應(yīng)用[J]. 工程力學(xué), 2014, 31(11): 39-45

Bai Junqiang, Wang Chen, Zhang Yang. Application of a Turbulence Model Based on Von-Karmen Length Scale in Steady and Unsteady Flow Simulation[J]. Engineering Mechanics, 2014, 31(11): 39-45 (in Chinese)

[13] Menter F R. A Scale-Adaptive Simiclation Modelusing Two-Equation Medels[R]. AIAA-2005-1095

[14] Menter F R. Two-Equation Eddy-Viscosity Turbulence Models for Egineering Applications[J]. AIAA Journal, 1994, 32(8): 1598-1605

[15] Sirovich L. Turbulence and the Dynamics of Coherent Structures[J]. Quarterly of Applied Mathematics, 1987, 45(3): 561-571

[16] Simpson T W, Mauery T M, Korte J J, et al. Kriging Models for Global Approximation in Simulation-Based Multidisciplinary Design Optimization[J]. AIAA Journal, 2001, 39(12): 2233-2241

[17] Li J, Huang S, Jiang S, et al. Unsteady Viscous Flow Simulations by a Fully Implicit Method with Deforming Grid[R]. AIAA-2005-1221