可變周期諧波平衡法求解周期性非定常渦脫落問(wèn)題*

柴振霞 劉偉 楊小亮 周云龍

(國(guó)防科技大學(xué)空天科學(xué)學(xué)院,長(zhǎng)沙 410073)

1 引 言

周期性非定常流動(dòng)是流體力學(xué)領(lǐng)域科學(xué)研究與實(shí)際工程應(yīng)用中常見(jiàn)的流動(dòng)現(xiàn)象,根據(jù)非定常性產(chǎn)生原因可以將其分為兩類(lèi)：第一類(lèi)是由物體周期性運(yùn)動(dòng)引起并強(qiáng)制驅(qū)動(dòng)的,這類(lèi)流動(dòng)都有一個(gè)明顯的特征頻率,物理量的波動(dòng)周期已知,如渦輪機(jī)械內(nèi)部轉(zhuǎn)子和靜子相對(duì)運(yùn)動(dòng)引起的非定常流動(dòng)、直升機(jī)的旋翼繞流、翼形的俯仰振蕩等;第二類(lèi)是由流動(dòng)本身的不穩(wěn)定性導(dǎo)致的,這類(lèi)流動(dòng)的波動(dòng)周期事先未知,如圓柱繞流中的渦脫落現(xiàn)象及其他包括分離和自由剪切層的流動(dòng)[1].

通常采用時(shí)域方法(time domain method,TDM)數(shù)值求解周期性非定常流動(dòng)問(wèn)題,即在時(shí)間方向上順序求解流動(dòng)控制方程,直到流動(dòng)達(dá)到周期性穩(wěn)定狀態(tài).時(shí)域計(jì)算方法思想簡(jiǎn)單,可用于各種非定常流動(dòng)問(wèn)題.但是,時(shí)域計(jì)算方法沒(méi)有利用流動(dòng)的周期特性,從初始條件迭代獲得穩(wěn)定周期解之前需要經(jīng)歷很長(zhǎng)階段的瞬態(tài)解的計(jì)算,而這部分的計(jì)算時(shí)間通常比真實(shí)波動(dòng)周期長(zhǎng)很多[2],因此其計(jì)算效率較低.同時(shí),時(shí)域計(jì)算通常使用較小的物理時(shí)間步長(zhǎng)以滿足穩(wěn)定性條件,從而導(dǎo)致非定常流動(dòng)計(jì)算的時(shí)間較長(zhǎng).

對(duì)于周期性流動(dòng),人們更關(guān)注的是流動(dòng)達(dá)到周期性的穩(wěn)定狀態(tài).近二十年來(lái),許多學(xué)者利用流動(dòng)變化的周期性特點(diǎn)發(fā)展了各種不同的頻域方法,如線性頻域方法、非線性諧波法、諧波平衡法、非線性頻域法和時(shí)間譜方法等.早期的線性頻域方法[3,4]把流動(dòng)變量分解為時(shí)均值和周期性小擾動(dòng)量的和,將流動(dòng)控制方程分解為不耦合的時(shí)均方程和擾動(dòng)方程,并基于小擾動(dòng)線性化假設(shè)對(duì)擾動(dòng)方程進(jìn)行線性化處理,通過(guò)求解線性化方程來(lái)模擬非定常擾動(dòng).該方法求解速度較快,但不適用于存在較強(qiáng)非線性的流動(dòng)問(wèn)題.為了能夠模擬非定常流場(chǎng)中的非線性效應(yīng),Ning和He[5]提出了非線性諧波法,該方法將時(shí)均方程和擾動(dòng)方程耦合求解.雖然這種方法能夠模擬一些變換復(fù)雜的非線性非定常流動(dòng),但仍然不能夠求解非簡(jiǎn)諧周期性流動(dòng).隨后,Hall等[6]提出了基于傅里葉展開(kāi)的諧波平衡法(harmonic balance method,HBM),該方法包含更多的非線性相互作用項(xiàng).由于非線性項(xiàng)在頻域內(nèi)的表示比較復(fù)雜,HBM不直接求解變量傅里葉級(jí)數(shù)中的諧波系數(shù),通過(guò)引入離散傅里葉變換,將非定常流動(dòng)求解過(guò)程轉(zhuǎn)換為一個(gè)周期內(nèi)幾個(gè)等時(shí)間間隔的瞬時(shí)流動(dòng)耦合求解.計(jì)算采用的時(shí)刻點(diǎn)越多,HBM計(jì)算包含的諧波分量越多,計(jì)算精度越高,因此該方法可以獲得較高的時(shí)間離散精度.Hall等[6]采用HBM成功模擬了單級(jí)、多級(jí)渦輪機(jī)葉片繞流問(wèn)題[7],計(jì)算表明即使是強(qiáng)非線性流動(dòng)HBM也能在較少的諧波數(shù)下達(dá)到工程要求的精度.McMullen等[8]發(fā)展了非線性頻域方法,該方法與HBM等價(jià),但在每一步迭代中都需要進(jìn)行快速傅里葉變換.Gopinath和Jameson[9]基于傅里葉變換進(jìn)一步提出了完全在時(shí)域求解流動(dòng)控制方程的時(shí)間譜方法,該方法在本質(zhì)上與HBM是相同的.時(shí)間譜方法被成功應(yīng)用到二維和三維非定常流動(dòng),如翼形俯仰振蕩和圓柱渦脫落問(wèn)題.

與TDM相比,HBM能夠在保證精度的同時(shí)大幅度地減少計(jì)算時(shí)間,在周期性的內(nèi)流和外流非定常問(wèn)題中得到了廣泛應(yīng)用.在內(nèi)流計(jì)算中,HBM被廣泛應(yīng)用于計(jì)算葉輪機(jī)械內(nèi)部的非定常流動(dòng)[10-17],如振蕩葉珊[6]、壓氣機(jī)轉(zhuǎn)靜干涉[13,14,17]、葉珊氣動(dòng)彈性顫振及極限環(huán)振動(dòng)[11,12]等問(wèn)題,與TDM相比,計(jì)算效率可實(shí)現(xiàn)1—2個(gè)數(shù)量級(jí)的提升.在外流計(jì)算中,Thomas等[18,19]和Guillaume等[20]對(duì)振蕩翼形和機(jī)翼的非定常計(jì)算表明,HBM計(jì)算準(zhǔn)確性高,可以模擬流動(dòng)中存在的激波等強(qiáng)非線性作用.Ekici等[21]采用HBM數(shù)值模擬了直升機(jī)旋翼的非定常繞流問(wèn)題.HBM在飛行器動(dòng)導(dǎo)數(shù)快速預(yù)測(cè)方面的應(yīng)用研究也取得了較好的效果[22-27].

由于HBM是在傅里葉級(jí)數(shù)展開(kāi)的基礎(chǔ)上構(gòu)建的,需要事先知道流動(dòng)穩(wěn)定周期解變化的頻率,因此HBM特別適用于已知頻率的第一類(lèi)周期性非定常流動(dòng)問(wèn)題.對(duì)于非定常周期事先不能確定的流動(dòng)問(wèn)題,如鈍體繞流中的渦脫落問(wèn)題,HBM計(jì)算可以采用實(shí)驗(yàn)測(cè)得的頻率或者TDM計(jì)算得到的頻率作為基準(zhǔn)頻率[28,29].由于鈍體繞流的非定常流動(dòng)是自發(fā)形成的,達(dá)到穩(wěn)定后的流動(dòng)變化周期與幾何布置和流動(dòng)參數(shù)有關(guān).因此,采用HBM計(jì)算時(shí)應(yīng)將頻率作為變量處理,以適用于不同幾何外形和來(lái)流條件.

為了求解這類(lèi)周期不可預(yù)知的非定常問(wèn)題,McMullen等[30]提出了一種基于殘差梯度的可變周期計(jì)算方法(gradient based variable time period method,GBVTP),并與非線性頻域方法相結(jié)合,實(shí)現(xiàn)了流動(dòng)變量和周期的同時(shí)迭代求解.McMullen等[30]采用該方法成功模擬了低雷諾數(shù)定姿態(tài)圓柱繞流中的渦脫落周期.Gopinath 和 Jameson[31]隨后提出了基于殘差梯度的可變周期時(shí)間譜方法,成功模擬了低雷諾數(shù)定姿態(tài)圓柱繞流和 NACA0012翼型大攻角繞流問(wèn)題.Spiker等[32]提出了一種基于相位差的方法確定圓柱繞流渦脫落頻率,并模擬了強(qiáng)迫振動(dòng)圓柱繞流下的渦脫落問(wèn)題.Mosahebi和Nadarajah[2,33]將GBVTP應(yīng)用到自適應(yīng)HBM.Yao和Jaiman[34]以及Yao和Marques[35]采用可變周期HBM數(shù)值模擬了低雷諾數(shù)圓柱渦致振蕩問(wèn)題以及跨聲速極限環(huán)振動(dòng)問(wèn)題,成功模擬了極限環(huán)振動(dòng)的振幅和頻率.國(guó)內(nèi)關(guān)于可變周期頻域方法的研究相對(duì)較少,比較典型的是2009年張煒和席光[36]采用基于殘差梯度的可變周期時(shí)間譜方法求解了二維不可壓方柱繞流問(wèn)題.總體來(lái)看,可變周期HBM的研究還比較少,需要開(kāi)展深入研究.

本文基于上述文獻(xiàn),進(jìn)一步開(kāi)展可變周期HBM在周期性非定常渦脫落問(wèn)題中的應(yīng)用研究.以二維圓柱和方柱繞流為例,詳細(xì)考察了該方法的計(jì)算精度和效率,并分析了不同參數(shù)對(duì)計(jì)算的影響.在GBVTP中,采用不同尋優(yōu)方法搜索周期T,對(duì)比研究了不同優(yōu)化方法的計(jì)算性能.

2 計(jì)算方法

2.1 流動(dòng)控制方程

以RANS方程為控制方程,其有限體積離散形式可以寫(xiě)為

其中J為雅克比行列式,計(jì)算網(wǎng)格單元體積V=J-1;Q=(ρ,ρu,ρv,ρw,ρe)T為流動(dòng)守恒變量;R(Q)為無(wú)黏通量E,F,G和黏性通量Ev,Fv,Gv空間離散產(chǎn)生的殘差向量,

本文無(wú)黏通量的空間離散采用AUSMPW+格式,黏性通量空間離散采用格林公式.

2.2 時(shí)域計(jì)算方法

非定常時(shí)域計(jì)算采用Jameson[37]提出的含雙時(shí)間步的LU-SGS隱式時(shí)間格式對(duì)控制方程進(jìn)行時(shí)間離散,推進(jìn)求解如下方程：

其中Δt為真實(shí)物理時(shí)間步長(zhǎng),Δτ為虛擬時(shí)間步長(zhǎng).

2.3 諧波平衡法

對(duì)于周期性流動(dòng),其流動(dòng)守恒變量Q和殘差R(Q)在時(shí)間上以頻率ω周期變化,可以表示為如下傅里葉級(jí)數(shù)的形式：

其中NH為諧波數(shù).

將方程(4)代入到方程(1)中,整理可得頻域形式諧波平衡方程.由于直接求解頻域方程十分困難,Hall等[6]將一個(gè)周期分為NT(NT=2NH+1)等份,利用離散傅里葉變換將頻域方程轉(zhuǎn)換回時(shí)域進(jìn)行求解.時(shí)域形式諧波平衡方程可以寫(xiě)為

詳細(xì)推導(dǎo)過(guò)程可以參考文獻(xiàn)[27].式中J-1ωDQ為諧波源項(xiàng),D是NT×NT系數(shù)矩陣,且

Q和R表示這NT個(gè)時(shí)刻的流動(dòng)變量和殘差向量,則有

式中ΔT=T/NT.

引入虛擬時(shí)間導(dǎo)數(shù)項(xiàng)來(lái)推進(jìn)求解方程(5)：

本文采用隱式LU-SGS方法求解方程(8),其中諧波源項(xiàng)顯式計(jì)算,推進(jìn)求解如下方程組：

其中Qi表示第ti=t0+iΔt時(shí)刻的守恒變量,

2.4 可變周期計(jì)算方法

時(shí)間頻率ω=2π/T,諧波平衡方程(8)可寫(xiě)為

鈍體繞流中的非定常性是自發(fā)形成的,在不同參數(shù)下,非定常流動(dòng)的周期未知.在非定常流動(dòng)的每一時(shí)刻,當(dāng)且僅當(dāng)T等于正確的周期T*時(shí)才能使控制方程組(10)相容,同時(shí)正確地反映流動(dòng)現(xiàn)象.對(duì)于時(shí)間周期不可預(yù)知的這類(lèi)周期性流動(dòng),本文采用基于殘差梯度的方法從初始猜測(cè)值開(kāi)始迭代求解正確的時(shí)間周期.

定義殘差

當(dāng)T=T*時(shí),殘差RHB=0,否則RHB≠0.因此,給定不同的周期T,HBM計(jì)算的殘差將收斂到不同精度水平.T越接近真實(shí)值,殘差越接近于零.可變周期HBM的目的是求解真實(shí)的周期T*使RHB→ 0,這是無(wú)約束的最優(yōu)化問(wèn)題,通常采用最速下降法(steepest descent method,SDM)進(jìn)行求解.

定義變量

對(duì)于殘差變量L,不同研究采用的定義不同,如文獻(xiàn)[36]基于y方向動(dòng)量方程的殘差構(gòu)造變量L.數(shù)值計(jì)算表明,不同的L定義對(duì)收斂速度影響較大,當(dāng)其量級(jí)較小時(shí),收斂較慢,計(jì)算時(shí)間長(zhǎng).計(jì)算表明,本文采用所有殘差分量平方之和的方式總體收斂速度較快.后續(xù)研究可以圍繞不同殘差分量的歸一化或者加權(quán)來(lái)構(gòu)造殘差變量L.

建立殘差梯度

利用殘差梯度更新時(shí)間周期：

選擇合適的初值T0和步長(zhǎng)λ,通過(guò)求解方程(10)和方程(14)流動(dòng)變量和周期T同時(shí)迭代更新并達(dá)到收斂.

周期搜索過(guò)程中,采用殘差的相對(duì)值判斷收斂,定義

式中RHS為所有時(shí)刻所有網(wǎng)格點(diǎn)所有殘差分量絕對(duì)值之和對(duì)總網(wǎng)格點(diǎn)的平均,RHS0為初始迭代的殘差.對(duì)于本文算例,ε≤10-3能夠滿足精度要求,因此計(jì)算采用ε=10—3.

3 數(shù)值方法驗(yàn)證

NACA0012翼型繞1/4弦點(diǎn)做俯仰振蕩是一個(gè)典型的周期性非定常問(wèn)題,經(jīng)常作為數(shù)值方法的驗(yàn)證算例.翼型的正弦俯仰振蕩運(yùn)動(dòng)采用如下攻角隨時(shí)間的變換規(guī)律：

式中α0為平均攻角,αm為振幅,為減縮頻率,為俯仰振動(dòng)的頻率,為來(lái)流速度,無(wú)量綱化長(zhǎng)度取翼型的弦長(zhǎng).本文計(jì)算選取AGARD實(shí)驗(yàn)中的CT5算例[38],主要的計(jì)算條件見(jiàn)表1.

表1 NACA0012翼型俯仰振蕩AGARD CT5算例計(jì)算條件Table 1.Computational conditions of the AGARD CT5 test case for the NACA0012 airfoil.

3.1 網(wǎng)格劃分

圖1是計(jì)算網(wǎng)格情況,采用O型拓?fù)浣Y(jié)構(gòu),網(wǎng)格數(shù)量為84×287 (法向×周向).為了保證遠(yuǎn)場(chǎng)邊界對(duì)計(jì)算的影響較小,遠(yuǎn)場(chǎng)邊界設(shè)在20倍弦長(zhǎng)以外.

圖1 NACA0012翼型計(jì)算網(wǎng)格Fig.1.Mesh for the NACA0012 airfoil.

3.2 計(jì)算結(jié)果及分析

采用HBM計(jì)算時(shí),諧波數(shù)NH=1,2,3,4,即分別將一個(gè)周期3,5,7,9等分,計(jì)算得到各個(gè)時(shí)刻點(diǎn)的瞬時(shí)流場(chǎng).采用TDM計(jì)算時(shí),無(wú)量綱時(shí)間步長(zhǎng)為0.01,子迭代收斂的判據(jù)為殘差下降兩個(gè)量級(jí).一個(gè)周期內(nèi)NACA0012翼型的升力系數(shù)和力矩系數(shù)隨攻角的變化如圖2所示.將HBM與TDM的計(jì)算結(jié)果與AGARD的實(shí)驗(yàn)數(shù)據(jù)[38]以及Batina[39]的計(jì)算結(jié)果進(jìn)行了比較.對(duì)于線性特征明顯的升力系數(shù)曲線,HBM在NH=1的計(jì)算結(jié)果與TDM計(jì)算結(jié)果符合良好,且與Batina[39]的計(jì)算結(jié)果一致,但計(jì)算結(jié)果比實(shí)驗(yàn)結(jié)果偏小.與升力系數(shù)相比,俯仰力矩系數(shù)具有較強(qiáng)的非線性,HBM計(jì)算至少需要3個(gè)諧波數(shù)才能模擬力矩系數(shù)隨攻角的變化規(guī)律.HBM重建的力矩系數(shù)曲線和時(shí)域計(jì)算結(jié)果及實(shí)驗(yàn)值能夠較好地符合,驗(yàn)證了HBM對(duì)周期性非定常運(yùn)動(dòng)的模擬能力.

圖2 NACA0012翼型的(a)升力系數(shù)和(b)俯仰力矩系數(shù)遲滯曲線Fig.2.(a)Lift and (b)pitching moment coefficients dynamic dependence of NACA0012 airfoil.

圖3展示了瞬時(shí)壓力系數(shù)沿翼型表面的分布規(guī)律.數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)[38]基本符合.HBM計(jì)算保留前3階諧波分量能夠較為準(zhǔn)確地模擬翼型表面壓力系數(shù)的變化規(guī)律,即使在激波出現(xiàn)的位置也具有較高的精度.

采用HBM計(jì)算時(shí),俯仰力矩系數(shù)收斂速度幾乎不受諧波數(shù)的影響,如圖4所示.HBM計(jì)算迭代3000步時(shí),動(dòng)態(tài)氣動(dòng)力矩系數(shù)達(dá)到收斂狀態(tài).為了定量分析HBM的計(jì)算效率,定義加速比圖5給出了加速度比隨諧波數(shù)的變化.其中,為了保證力矩系數(shù)的收斂性,TDM模擬兩個(gè)周期的俯仰運(yùn)動(dòng).HBM計(jì)算所需CPU時(shí)間按3000步計(jì)算.由圖5可知,保留1階諧波分量時(shí),HBM的計(jì)算效率比TDM高一個(gè)量級(jí),隨著諧波數(shù)的增加,HBM的計(jì)算效率下降.諧波數(shù)等于3時(shí),HBM計(jì)算比TDM大約快4倍.

圖3 NACA0012翼型俯仰振蕩過(guò)程中的瞬時(shí)壓力系數(shù)分布 (a)攻角減小過(guò)程中α=—2.41°;(b)攻角增大過(guò)程中α=—2.00°Fig.3.Instantaneous pressure coefficient distribution compared to experimental data of NACA0012 airfoil：(a)α=—2.41° for decreasing angle;(b)α=—2.00° for increasing angle.

圖4 HBM取不同諧波數(shù)時(shí)俯仰力矩系數(shù)收斂曲線 (a)NH=1;(b)NH=3Fig.4.Pitching moment coefficient convergence history for the HBM with respect to the number of harmonics：(a)NH=1;(b)NH=3.

圖5 CPU時(shí)間加速比隨諧波數(shù)的變化Fig.5.CPU time speedup of the HBM with respect to the TDM.

4 非定常圓柱繞流中的渦脫落問(wèn)題

非定常圓柱繞流的流態(tài)與雷諾數(shù)密切相關(guān),當(dāng)雷諾數(shù)47 ＜Re＜ 194時(shí),圓柱底部會(huì)出現(xiàn)二維的周期性渦脫落,稱(chēng)之為卡門(mén)渦街.本文模擬低雷諾數(shù)下的二維非定常圓柱繞流,來(lái)流馬赫數(shù)為0.2,將模擬結(jié)果與Williamson的實(shí)驗(yàn)數(shù)據(jù)[43,44]及McMullen等[30]和Spiker等[32]的數(shù)值計(jì)算結(jié)果進(jìn)行比較.采用O型計(jì)算網(wǎng)格,如圖6所示,計(jì)算網(wǎng)格大小為105×61 (流向×法向),壁面第一層網(wǎng)格法向尺度為2.0×10—4D,D為圓柱直徑.

4.1 非定常時(shí)域計(jì)算

首先采用本課題組開(kāi)發(fā)的非定常時(shí)域計(jì)算軟件ADCRP,對(duì)來(lái)流Re=180條件下的圓柱繞流進(jìn)行數(shù)值模擬,計(jì)算得到升力系數(shù)CL和阻力系數(shù)隨時(shí)間的變化過(guò)程見(jiàn)圖7.計(jì)算得到的平均阻力系數(shù)CD0=1.3457,最大升力系數(shù)CLmax=0.6347,根據(jù)升力系數(shù)振動(dòng)頻率得到渦脫落頻率St=0.185,對(duì)應(yīng)無(wú)量綱周期T=5.41.計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合良好,如表2所列.

圖7 升、阻力系數(shù)收斂曲線Fig.7.Time history of lift coefficientCLand dragCD.

表2 時(shí)域計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比Table 2.Time-averaged coefficient and Strouhal number compared with experiment data.

4.2 可變周期諧波平衡法計(jì)算

4.2.1 諧波數(shù)收斂性分析

首先,采用可變周期HBM對(duì)來(lái)流Re=180條件下的圓柱繞流進(jìn)行數(shù)值模擬,考察計(jì)算所需諧波數(shù).初始周期T0=4,步長(zhǎng)λ=1.

圖8為不同諧波數(shù)下計(jì)算得到的周期T收斂曲線,圖9為諧波數(shù)取3時(shí)不同時(shí)刻升力系數(shù)收斂曲線.由圖8和圖9可知,非定常渦脫落周期和氣動(dòng)力系數(shù)同時(shí)迭代更新并達(dá)到收斂.圖10和圖11分別為不同諧波數(shù)下重建得到的升力系數(shù)和阻力系數(shù)在一個(gè)周期內(nèi)的變化曲線.表3列出了采用不同諧波數(shù)計(jì)算得到的渦脫落頻率和平均阻力系數(shù).從計(jì)算結(jié)果可以看出,隨著諧波數(shù)的增加,氣動(dòng)力及渦脫落周期逐漸收斂.諧波數(shù)NH=3計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及TDM計(jì)算結(jié)果相符合,諧波數(shù)再增加時(shí)計(jì)算精度沒(méi)有明顯變化,說(shuō)明3個(gè)諧波數(shù)能夠滿足計(jì)算精度要求.

圖8 不同諧波數(shù)下的周期T收斂曲線Fig.8.Convergence from initial guess to exact time period with varying number of harmonics.

圖12為流動(dòng)變化周期內(nèi)3個(gè)等間隔時(shí)刻的流線圖,清晰地表現(xiàn)出圓柱后緣交替產(chǎn)生的渦脫落現(xiàn)象.圖13為HBM還原得到的升力系數(shù)最小時(shí)刻的非定常瞬態(tài)熵流場(chǎng)與時(shí)域計(jì)算結(jié)果的對(duì)比.采用3階諧波計(jì)算,就可以很好地還原圓柱非定常渦脫落流場(chǎng)的多層次渦結(jié)構(gòu).

圖9 升力系數(shù)收斂曲線(NH=3)Fig.9.Time history of lift coefficientCLwithNH=3.

圖10 升力系數(shù)隨時(shí)間的變化Fig.10.Variation ofCLover one period.

圖11 阻力系數(shù)隨時(shí)間的變化Fig.11.Variation ofCDover one period.

取3階諧波對(duì)雷諾數(shù) 60≤Re≤180 的圓柱繞流進(jìn)行數(shù)值模擬,并將計(jì)算得到的斯特勞哈爾數(shù)St和CD0與Williamson[43,44]和Henderson[40]的實(shí)驗(yàn)數(shù)據(jù)及McMullen等[30]和Spiker等[32]的數(shù)值計(jì)算結(jié)果相比較.圖14為表征非定常運(yùn)動(dòng)周期特征的St隨雷諾數(shù)的變化,圖15為平均阻力系數(shù)CD0隨雷諾數(shù)的變化.本文可變周期HBM在3階諧波數(shù)下的計(jì)算結(jié)果與實(shí)驗(yàn)值符合良好,誤差在4%以?xún)?nèi).對(duì)于St,數(shù)值計(jì)算結(jié)果都低于實(shí)驗(yàn)值,本文計(jì)算結(jié)果與其他數(shù)值計(jì)算結(jié)果一致.對(duì)于平均阻力系數(shù),可變周期HBM的計(jì)算結(jié)果比McMullen等[30]采用可變周期非線性頻域方法計(jì)算的結(jié)果更接近于實(shí)驗(yàn)值.

表3 不同諧波數(shù)下的計(jì)算結(jié)果Table 3.Strouhal number and time-averaged coefficient computed by different number of harmonics.

圖12 Re=180,NH=3條件下不同時(shí)刻的流線圖 (a)t=T/3;(b)t=2T/3;(c)t=TFig.12.Streamlines at various time instances over one period (Re=180,NH=3)：(a)t=T/3;(b)t=2T/3;(c)t=T.

圖13 熵等值線圖(CL最小時(shí)刻)(a)TDM計(jì)算結(jié)果;(b)HBM計(jì)算結(jié)果（NH=3）Fig.13.Comparison of instantaneous entropy contours：(a)TDM results;(b)HBM results (NH=3).

圖14 Strouhal數(shù)隨Re的變化Fig.14.Strouhal number as a function of Reynolds number.

4.2.2 迭代步長(zhǎng)λ的影響分析

圖15 平均阻力系數(shù)隨Re的變化Fig.15.Mean coefficient of drag versus Reynolds number.

圖16為Re=180,NH=3情況下,取不同步長(zhǎng)λ計(jì)算得到的周期T.λ對(duì)最終結(jié)果影響較小,非定常渦脫落周期都收斂到同一值.λ=1時(shí)周期T收斂最慢,步長(zhǎng)增大收斂速度有所加快,但λ＞10之后收斂曲線沒(méi)有明顯變化.

圖16 不同步長(zhǎng)λ下的周期T收斂曲線Fig.16.Time period convergence computed with three different step sizesλ.

4.2.3 初始周期T0對(duì)計(jì)算的影響分析

圖17為Re=180,NH=3情況下,取不同初始T0計(jì)算得到的周期T.當(dāng)T0＜ 5.8時(shí),取不同初值T0,渦脫落周期最終都收斂到真實(shí)周期T*=5.389;當(dāng) 5.8≤T0≤13 時(shí),取不同初值T0,最終結(jié)果都收斂到T=11.43,與T*近似成兩倍的關(guān)系,升力系數(shù)在 0≤t≤11.43 內(nèi)也確實(shí)變化兩個(gè)周期,如圖18所示;當(dāng)T0再增大到17和20的時(shí)候,渦脫落周期最終都收斂到T=16.834,與T*近似成三倍的關(guān)系.對(duì)于周期性非定常問(wèn)題,如果T*是流動(dòng)變化周期,那么2T*和3T*也是流動(dòng)變化的周期.對(duì)于可變周期HBM,在周期T搜索過(guò)程中有可能收斂到nT*(n=2,3,···).值得注意的是,實(shí)際計(jì)算結(jié)果并不是準(zhǔn)確的倍數(shù)關(guān)系,這可能與諧波數(shù)有關(guān),也可能是計(jì)算誤差導(dǎo)致的,此現(xiàn)象值得進(jìn)一步深入研究.同時(shí),有必要發(fā)展自動(dòng)搜索最小周期的計(jì)算方法.

4.2.4 計(jì)算效率分析

圖19比較了采用HBM與TDM計(jì)算得到的St隨雷諾數(shù)的變化曲線.當(dāng)時(shí)域計(jì)算選擇時(shí)間步長(zhǎng)為Δt=0.01時(shí),在雷諾數(shù) 8 0≤Re≤180 之間兩種計(jì)算方法的結(jié)果符合良好,雷諾數(shù)低于80時(shí),HBM的計(jì)算結(jié)果比TDM結(jié)果更接近實(shí)驗(yàn)值.但當(dāng)TDM計(jì)算采用較小步長(zhǎng)Δt=0.001時(shí),較低雷諾數(shù)下的TDM計(jì)算結(jié)果與HBM計(jì)算保留前3階諧波分量的結(jié)果基本重合,不同的是此時(shí)的TDM付出的計(jì)算代價(jià)更大.圖20為雷諾數(shù)Re=180和Re=60的加速比.由圖20可知,來(lái)流雷諾數(shù)越低,HBM相對(duì)TDM的計(jì)算效率優(yōu)勢(shì)越明顯.

圖17 不同步長(zhǎng)T0下計(jì)算的周期T收斂曲線Fig.17.Time period convergence with various starting guessesT0.

圖18T=11.43 時(shí)重建的升力系數(shù)曲線Fig.18.Variation ofCLover one period with converged time periodT=11.43.

圖19 HBM計(jì)算的St與TDM計(jì)算結(jié)果的對(duì)比Fig.19.Comparison of the HBMStdata results with TDM results.

圖20 不同雷諾下的加速比Fig.20.CPU time speedup of various Reynolds number.

4.3 基于固定周期諧波平衡法的渦脫落周期辨識(shí)

對(duì)于固定周期HBM計(jì)算,當(dāng)給定的周期T不等于真實(shí)的物理周期T*時(shí),殘差將下降到一定水平,并且收斂后的殘差和氣動(dòng)力會(huì)隨著虛擬迭代步的增加呈現(xiàn)周期性變化,如圖21所示,其中ρ為密度;u,v分別為x方向和y方向的速度;e為單位質(zhì)量總能.當(dāng)Re=180時(shí),采用可變周期HBM計(jì)算的圓柱繞流渦脫落周期為T(mén)*=5.389.固定周期HBM給定T=4時(shí),計(jì)算得到的升力系數(shù)和殘差都呈現(xiàn)周期性變化的趨勢(shì),相鄰迭代步重建的升力系數(shù)曲線之間存在相位差,如圖22所示.對(duì)于給定的T,當(dāng)升力系數(shù)收斂后每一步虛擬迭代產(chǎn)生的相位差是常值.而且這個(gè)相位差與給定的周期T呈近似線性的關(guān)系,當(dāng)且僅當(dāng)給定的T為真實(shí)的周期T*時(shí)相位差為0,HBM計(jì)算得到各個(gè)時(shí)刻的升力系數(shù)曲線不再周期性振蕩,而是收斂到常值,如圖23所示.Spiker等[32]提出的相位差方法正是基于這一點(diǎn)來(lái)辨識(shí)渦脫落的頻率.圖24給出了固定周期HBM計(jì)算時(shí),相位差隨周期的變化.由圖24可知,在一定范圍內(nèi),相位差的確與T成線性關(guān)系,采用固定周期HBM辨識(shí)渦脫落周期時(shí),只需計(jì)算兩個(gè)不同周期對(duì)應(yīng)的相位差就可以通過(guò)線性關(guān)系推算相位差等于0時(shí)對(duì)應(yīng)的T*.值得注意的是,采用3階諧波進(jìn)行計(jì)算時(shí),當(dāng)T=11.43時(shí)相位差也等于0,在該值附近也存在線性關(guān)系,說(shuō)明11.43也是旋渦脫落的周期,這與之前可變周期HBM的計(jì)算結(jié)果一致.

圖21 升力系數(shù)和t=T時(shí)刻的殘差收斂曲線(Re=180,T=4,NH=3)(a)升力系數(shù);(b)殘差Fig.21.Time history of lift coefficientCLat various time instances over one period and residual att=T(Re=180,T=4,NH=3)：(a)Lift coefficient;(b)residual.

圖22 不同迭代步重建的升力系數(shù)隨時(shí)間的變化(Re=180,T=4,NH=3)(a)整體;(b)局部Fig.22.Variation ofCLover one period at different iterations：(a)Overall;(b)local.

圖23 T=5.389時(shí)各個(gè)時(shí)刻升力系數(shù)收斂曲線(NH=3)Fig.23.Time history of lift coefficientCLat various time instances over one period withT=5.389 (NH=3).

GBVTP的實(shí)現(xiàn)是基于殘差在T=T*時(shí)最小的物理事實(shí),采用SDM沿著殘差下降的方向進(jìn)行搜索.圖25為采用固定周期HBM計(jì)算的殘差隨周期T的變化.由圖25可知,T的精度決定了HBM計(jì)算的精度.當(dāng)給定的周期遠(yuǎn)離T*時(shí),殘差下降越來(lái)越困難,當(dāng)T=T*時(shí)殘差迅速下降.因?yàn)?當(dāng)且僅當(dāng)給定的T為真實(shí)的渦脫落周期才能滿足非定常流動(dòng)控制方程,使得殘差接近于0.因此,采用3個(gè)諧波數(shù)進(jìn)行計(jì)算時(shí),T=5.389和11.43都是旋渦脫落的周期,這與可變周期計(jì)算的結(jié)果也相一致.

圖24 相位差隨周期T的變化(Re=180,NH=3)Fig.24.Change in phase of unsteady lift versus time period forRe=180 (NH=3).

圖25 殘差隨周期T的變化(Re=180,NH=3)Fig.25.HBM solution residual versus time period forRe=180 (NH=3).

5 多種尋優(yōu)策略在可變周期計(jì)算方法中的比較

GBVTP的核心思想是尋找T*使得殘差趨于0,這是單變量單目標(biāo)函數(shù)的無(wú)約束最優(yōu)化問(wèn)題.求解這類(lèi)問(wèn)題的最優(yōu)化方法大致分為兩類(lèi)[45]：一類(lèi)計(jì)算用到函數(shù)的導(dǎo)數(shù);一類(lèi)只用目標(biāo)函數(shù)不計(jì)算導(dǎo)數(shù),為直接方法.這里考察了用到函數(shù)導(dǎo)數(shù)的另外兩種方法：牛頓法和Fletcher-Reeves共軛梯度法(簡(jiǎn)稱(chēng)FR法),并和前面計(jì)算采用的SDM進(jìn)行了比較.

5.1 優(yōu)化方法

考慮無(wú)約束問(wèn)題：

其中RHB表示諧波平衡方程的殘差向量.這個(gè)問(wèn)題是求解T*使得L(T*)→ 0.

采用SDM計(jì)算時(shí),沿著L下降最快的方向(負(fù)梯度方向)尋找T*,迭代公式為

其中λ表示搜索步長(zhǎng).

采用牛頓法計(jì)算時(shí),沿著牛頓方向計(jì)算T*,迭代公式為

采用FR法計(jì)算時(shí),需要構(gòu)造一組共軛方向,然后沿著這組方向進(jìn)行搜索,求解T*：

不可忽略的是,在FR方法中,初始的搜索方向必須選擇最速下降方向,即：

5.2 計(jì)算結(jié)果及分析

分別采用以上方法對(duì)Re=180的圓柱繞流進(jìn)行數(shù)值模擬,對(duì)比他們的計(jì)算精度和效率.圖26給出了牛頓法及SDM的計(jì)算結(jié)果,可以看出,牛頓法模擬的渦脫落周期與SDM的結(jié)果一致.采用SDM時(shí),需要設(shè)置搜索步長(zhǎng)λ,步長(zhǎng)越大收斂越快.而采用牛頓法計(jì)算時(shí)不需要設(shè)置此參數(shù),其收斂曲線介于SDM采用λ=1和λ=100的收斂曲線之間,且接近于λ=100的計(jì)算結(jié)果.

與SDM類(lèi)似,共軛梯度法也需要設(shè)置搜索步長(zhǎng)λ.在本文算例中,取不同時(shí)間步長(zhǎng)計(jì)算時(shí),共軛梯度法計(jì)算得到的周期T收斂曲線基本重合,步長(zhǎng)對(duì)共軛梯度法計(jì)算結(jié)果的影響不大,如圖27所示.共軛梯度法的收斂速度與SDM采用λ=100時(shí)的相同.

圖26 采用牛頓法和SDM計(jì)算的周期T收斂曲線對(duì)比圖 (a)初始T0=4;(b)初始T0=5.41Fig.26.Convergence of shedding time period computed by Newton method and SDM：(a)T0=4;(b)T0=5.41.

可變周期HBM采用三種不同優(yōu)化方法計(jì)算得到的周期T收斂曲線對(duì)比如圖28所示.由圖28可知,三種方法的計(jì)算精度相同,最終計(jì)算得到的渦脫落周期T均為5.389.共軛梯度法和牛頓法的收斂速度與SDM采用最大搜索步長(zhǎng)時(shí)的收斂速度一致.由于牛頓法不需要設(shè)置搜索步長(zhǎng)等參數(shù),在實(shí)際工程計(jì)算中更有優(yōu)勢(shì).

圖27 采用FR法計(jì)算的周期T收斂曲線(a)及其與SDM計(jì)算結(jié)果的比較(b)Fig.27.Convergence of shedding time period computed by FR conjugate gradient method (a)and compared with the SDM results (b).

圖28 采用三種不同優(yōu)化方法計(jì)算得到的周期T收斂曲線圖Fig.28.Convergence of shedding time period computed by three different methods of optimization.

6 非定常方柱繞流中的渦脫落問(wèn)題

采用可變周期HBM數(shù)值模擬二維不可壓非定常方柱繞流.來(lái)流雷諾數(shù)為Re=100,采用H型計(jì)算網(wǎng)格(如圖29所示),網(wǎng)格大小為147×106.為了對(duì)比研究,時(shí)域計(jì)算方法的結(jié)果列于表4.計(jì)算結(jié)果表明,TDM計(jì)算采用時(shí)間步長(zhǎng)0.01能夠滿足計(jì)算精度要求.

表5列出了渦脫落頻率、平均阻力系數(shù)和加速比隨諧波數(shù)的變化.隨著諧波數(shù)的增加,計(jì)算結(jié)果逐漸收斂.對(duì)于渦脫落頻率,采用3個(gè)諧波的計(jì)算結(jié)果與4個(gè)諧波的相同,且與TDM計(jì)算及實(shí)驗(yàn)值相符合.與TDM相比,HBM計(jì)算保留3個(gè)諧波數(shù)時(shí)計(jì)算效率提高了將近18倍.諧波數(shù)越少,計(jì)算速度越快,但精度會(huì)有所損失.當(dāng)諧波數(shù)增加時(shí),需要減小虛擬時(shí)間步長(zhǎng)以滿足計(jì)算穩(wěn)定性的要求,因此收斂速度明顯減慢.一個(gè)周期內(nèi)升力系數(shù)對(duì)比如圖30所示.圖31為升力系數(shù)最小時(shí)刻對(duì)應(yīng)的熵等值線圖.從圖30和圖31可以看出,HBM采用3個(gè)諧波數(shù)就能夠重現(xiàn)時(shí)域結(jié)果,達(dá)到時(shí)域計(jì)算的精度.

圖29 二維方柱繞流計(jì)算網(wǎng)格Fig.29.Computational grid for rectangular in cross flow.

表4 時(shí)域計(jì)算結(jié)果Table 4.Time-averaged coefficient and Strouhal number computed by time-domain solver using different physical time steps.

表5 Re=100時(shí)不同諧波數(shù)下的計(jì)算結(jié)果對(duì)比Table 5.Convergency of frequency and time-averaged coefficient with speedup estimates.

圖30 升力系數(shù)隨時(shí)間的變化Fig.30.Comparison of lift coefficients of HBM and TDM at Re=100.

圖31 熵等值線圖(CL最小時(shí)刻)(a)TDM計(jì)算結(jié)果;(b)HBM計(jì)算結(jié)果(NH=3)Fig.31.Comparison of the instantaneous entropy contours：(a)TDM results;(b)HBM results (NH=3).

7 結(jié) 論

本文以RANS方程為控制方程,采用HBM對(duì)低雷諾數(shù)二維不可壓繞圓柱和方柱的周期性非定常流動(dòng)進(jìn)行了數(shù)值模擬.對(duì)于這類(lèi)流動(dòng)周期未知的非定常流動(dòng)問(wèn)題,采用基于殘差導(dǎo)數(shù)的GBVTP求解渦脫落周期.得到的以下主要結(jié)論.

1)可變周期HBM可以準(zhǔn)確地模擬周期性非定常渦脫落問(wèn)題,計(jì)算得到的Strouhal數(shù)和平均阻力系數(shù)與實(shí)驗(yàn)值及其他數(shù)值計(jì)算數(shù)據(jù)符合良好,與傳統(tǒng)的TDM相比該方法具有較高的計(jì)算效率.

2)周期搜索步長(zhǎng)λ對(duì)計(jì)算結(jié)果影響較小,非定常渦脫落周期都收斂到同一值;當(dāng)初值T0較大時(shí),最終計(jì)算得到的周期T可能會(huì)收斂到nT*(n=2,3,···),因此有必要發(fā)展自動(dòng)搜索最小周期的計(jì)算方法.

3)基于不同優(yōu)化策略的可變周期計(jì)算結(jié)果表明：不同優(yōu)化方法的計(jì)算精度相當(dāng);牛頓法沒(méi)有參數(shù)問(wèn)題,其收斂速度與共軛梯度法及SDM采用最大搜索步長(zhǎng)時(shí)的收斂速度一致.因此,牛頓法在工程計(jì)算中更有優(yōu)勢(shì).