柴振霞 劉偉 楊小亮 周云龍

(國防科技大學(xué)空天科學(xué)學(xué)院,長沙 410073)

1 引 言

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

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

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

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

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

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

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

2 計算方法

2.1 流動控制方程

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

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

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

2.2 時域計算方法

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

其中Δt為真實物理時間步長,Δτ為虛擬時間步長.

2.3 諧波平衡法

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

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

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

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

Q和R表示這NT個時刻的流動變量和殘差向量,則有

式中ΔT=T/NT.

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

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

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

2.4 可變周期計算方法

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

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

定義殘差

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

定義變量

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

建立殘差梯度

利用殘差梯度更新時間周期：

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

周期搜索過程中,采用殘差的相對值判斷收斂,定義

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

3 數(shù)值方法驗證

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

式中α0為平均攻角,αm為振幅,為減縮頻率,為俯仰振動的頻率,為來流速度,無量綱化長度取翼型的弦長.本文計算選取AGARD實驗中的CT5算例[38],主要的計算條件見表1.

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

3.1 網(wǎng)格劃分

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

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

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

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

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

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

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

圖3 NACA0012翼型俯仰振蕩過程中的瞬時壓力系數(shù)分布 (a)攻角減小過程中α=—2.41°;(b)攻角增大過程中α=—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ù)收斂曲線 (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ù)的變化Fig.5.CPU time speedup of the HBM with respect to the TDM.

4 非定常圓柱繞流中的渦脫落問題

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

4.1 非定常時域計算

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

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

表2 時域計算結(jié)果與實驗結(jié)果對比Table 2.Time-averaged coefficient and Strouhal number compared with experiment data.

4.2 可變周期諧波平衡法計算

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

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

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

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

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

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

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

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

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

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

圖12 Re=180,NH=3條件下不同時刻的流線圖 (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最小時刻)(a)TDM計算結(jié)果;(b)HBM計算結(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 迭代步長λ的影響分析

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

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

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

4.2.3 初始周期T0對計算的影響分析

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

4.2.4 計算效率分析

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

圖17 不同步長T0下計算的周期T收斂曲線Fig.17.Time period convergence with various starting guessesT0.

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

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

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

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

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

圖21 升力系數(shù)和t=T時刻的殘差收斂曲線(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ù)隨時間的變化(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ù)收斂曲線(NH=3)Fig.23.Time history of lift coefficientCLat various time instances over one period withT=5.389 (NH=3).

GBVTP的實現(xiàn)是基于殘差在T=T*時最小的物理事實,采用SDM沿著殘差下降的方向進(jìn)行搜索.圖25為采用固定周期HBM計算的殘差隨周期T的變化.由圖25可知,T的精度決定了HBM計算的精度.當(dāng)給定的周期遠(yuǎn)離T*時,殘差下降越來越困難,當(dāng)T=T*時殘差迅速下降.因為,當(dāng)且僅當(dāng)給定的T為真實的渦脫落周期才能滿足非定常流動控制方程,使得殘差接近于0.因此,采用3個諧波數(shù)進(jìn)行計算時,T=5.389和11.43都是旋渦脫落的周期,這與可變周期計算的結(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)策略在可變周期計算方法中的比較

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

5.1 優(yōu)化方法

考慮無約束問題：

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

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

其中λ表示搜索步長.

采用牛頓法計算時,沿著牛頓方向計算T*,迭代公式為

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

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

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

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

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

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

圖27 采用FR法計算的周期T收斂曲線(a)及其與SDM計算結(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)化方法計算得到的周期T收斂曲線圖Fig.28.Convergence of shedding time period computed by three different methods of optimization.

6 非定常方柱繞流中的渦脫落問題

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

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

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

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

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

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

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

7 結(jié) 論

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

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

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

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