諧波平衡法在低速非定常流模擬中的應用

杜鵬程 寧方飛

(北京航空航天大學 能源與動力工程學院，北京 100191)

諧波平衡法在低速非定常流模擬中的應用

杜鵬程 寧方飛

(北京航空航天大學 能源與動力工程學院，北京 100191)

諧波平衡法是一種有效的周期性非定常流的計算方法.采用基于可壓縮流的諧波平衡方程在計算低速不可壓流動時，會由于對流通量計算格式中的數(shù)值粘性污染，降低解的精度和收斂性.采用預處理技術，使得基于可壓縮流的諧波平衡方程可以直接用于低速周期性非定常流的計算中.選取典型的不可壓方腔驅(qū)動流和低雷諾數(shù)圓柱繞流為例進行了時間推進法和諧波平衡法的計算對比.計算結(jié)果表明預處理后的諧波平衡方程適合于低速流的計算，在諧波平衡法中采用較少階數(shù)的諧波計算就可以還原出幾乎準確的非定常流場.

周期性非定常流;諧波平衡法;低速預處理;方腔驅(qū)動流;圓柱繞流

周期性非定常流是一種在內(nèi)流、外流中常見的流動，如葉輪機的動靜干涉、震蕩葉柵、翼形的俯仰振蕩等.對于這類問題的求解通常采用傳統(tǒng)的時間推進法.基于時間推進法的非定常計算往往十分耗時，計算的大部分時間用于過渡到周期性狀態(tài)的中間計算過程.為了減少計算時間，傅里葉方法被引入到上述周期性非定常流計算中［1］.在這類基于傅里葉方法的頻域方法中，參考文獻［2］提出的諧波平衡法將頻域方程變換回時域求解，同時可以很方便的處理湍流模型問題，能在保證計算精度的同時，加快計算速度，從而得到了廣泛關注和應用［3］.

在低速不可壓流計算中，如果采用基于可壓縮流動方程的計算格式，則會由于對流通量計算格式中的數(shù)值粘性過大和數(shù)值剛性問題造成解的精度降低且收斂性差［4］.因此不能將基于可壓縮流的諧波平衡方程直接應用于低速周期性非定常流計算中.而時間預處理方法可有效地解決這一問題.該方法通過在可壓流控制方程的時間導數(shù)項左乘一個預處理矩陣，改變了對流通量系數(shù)矩陣的特征值，降低了其系數(shù)矩陣的條件數(shù)，改善了數(shù)值剛性;并相應修改其對流通量格式中的耗散項，從而可根據(jù)流動條件減小數(shù)值粘性.當然，預處理矩陣形式并不唯一，選定不同的預處理矩陣可推導出不同對流格式的預處理形式.參考文獻［5］對不同的預處理矩陣做了分析對比.

本文采用低速預處理，使得可以將基于可壓縮流的諧波平衡方程直接用于低速周期性非定常流的計算中.并用方腔驅(qū)動流和低雷諾數(shù)圓柱繞流這兩個基本算例進行了驗證.

1 流動控制方程

1.1 時間推進法

采用有限體積法離散的Navier-Stokes方程可以寫為

式中，Q 為守恒變量;Q=［ρ，ρu，ρv，ρw，ρE］T;V 為控制體的體積;R為殘差向量.對于定常計算，本文采用隱式算法LU-SGS求解，即是求解:

非定常計算采用雙時間步方法DTS(Dual Time Stepping)推進求解如下方程:

式中，Δt*為物理時間步長.

1.2 諧波平衡法

作為時間推進法的一種降階模型，諧波平衡法HB(Harmonic Balance method)可以計算時間周期性流動［2］.對于周期性流動，其解Q和殘差R(Q)同樣是時間的周期函數(shù)，頻率為ω，可以將其用有限階傅里葉級數(shù)近似表示為

將式(4)代入控制方程式(1)，利用正弦函數(shù)的正交性，對其進行諧波平衡，可得

式(5)為方程數(shù)和未知數(shù)都為NT=2NH+1個的方程組，寫成矩陣形式為

由于控制方程式(6)在頻域內(nèi)的形式很復雜，直接計算十分麻煩.文獻［2］巧妙的利用離散傅里葉變換將頻域方程(6)再轉(zhuǎn)化回時域中求解.取一個周期T=2π/ω內(nèi)均勻分布的NT個時刻的解，定義:

其中Δt=T/NT，利用離散傅里葉變換，可得

將式(8)帶入式(6)可得

其中D=E－1AE.對于式(9)的求解，引入虛擬時間推進求解:

可以看到，諧波平衡法相當于是同時求解一個周期內(nèi)NT個時刻的“定常”流場，并通過源項Shb將這NT個時刻的流場關聯(lián)起來.但該源項的存在會影響計算的穩(wěn)定性，使得計算可取的CFL數(shù)變小.本文諧波平衡方程的求解采用參考文獻［6］提出的隱式BJ-SSOR算法.

1.3 帶低速預處理的諧波平衡方程

諧波平衡法的控制方程(10)和控制方程(1)相比，其對流通量的求法不變.因此在計算低速流時采用基于壓縮流動方程所發(fā)展的對流格式，則會由于數(shù)值粘性過大和數(shù)值剛性問題造成解的精度低且收斂性差.本文采用預處理方法來解決這一問題.預處理的基本原理和預處理矩陣的具體推導參見文獻［5］，最終得到帶預處理的諧波平衡方程為

其中

2 算例1——方腔驅(qū)動流

2.1 定常計算

方腔驅(qū)動流通常用來進行低速不可壓流的計算驗證［8］.即是一個長度為D的方腔，上壁面以速度U做勻速運動，其余壁面靜止，由于流體的粘性，會在方腔內(nèi)形成一個穩(wěn)定的漩渦.本文取U=0.1m/s，保證雷諾數(shù) Re=100，流動的馬赫數(shù)約為 0.001.

在定常計算時，網(wǎng)格大小為均勻分布的41×41.不采用預處理時，由于數(shù)值粘性過大，使得方腔中的流體幾乎靜止.當選取預處理馬赫數(shù)Mr=0.001后，計算得到的方腔中間截面的速度分布和參考文獻［9］的計算結(jié)果吻合很好，如圖1所示.同樣，能準確計算得從右壁面－下壁面－左壁面的靜壓升系數(shù)分布，如圖2所示.

2.2 非定常計算

使上壁面做周期性往復運動:

本文取ω=1 rad/s.雙時間步非定常計算時一個周期內(nèi)采用60個物理時間步，每個物理時間步內(nèi)沿虛擬時間步推進求解30次，使其殘差下降1×10－3量級.沿時間推進求解3個周期左右，流動就已達到良好的周期性狀態(tài)，上壁面阻力系數(shù)隨時間的變化過程如圖3所示.

圖1 方腔中間截面的速度分布

圖2 沿壁面的靜壓升系數(shù)分布

圖3 時間推進法計算的阻力系數(shù)隨時間的變化

在諧波平衡法計算中分別采用了不同階數(shù)諧波，并還原得到t=0.5T時刻的流場，方腔內(nèi)的流線與時間推進法計算結(jié)果的對比如圖4所示.雖然上壁面做一階諧波非定常運動，但是由于流體控制方程的非線性，如果僅僅采用一階諧波，還原得到的流場誤差較大，該時刻渦的運動中心有明顯區(qū)別，位置更加靠近中間截面.隨著采取的諧波階數(shù)的提高，還原得到的流場更加準確.采用3階諧波，即計算7個時刻的流場，就可以很好的還原出該時刻的非定常流場.

計算得到的一個周期內(nèi)上壁面的阻力系數(shù)的變化如圖5所示.由于阻力系數(shù)的計算只涉及到靠近上壁面處的流場，采用1階諧波就可以得到較準確的阻力系數(shù)的分布.當采用的諧波階數(shù)N≥3后，再增加諧波的階數(shù)對于解的精度沒有明顯的提高，并和時間推進法計算的阻力系數(shù)的幅值存在較小的差別.這是因為在對流通量和粘性通量這些非線性項計算時，采用有限階的諧波截斷和傅里葉變換總是拋棄了一些高階項，存在一定的計算誤差.

圖4 t=0.5T時刻方腔內(nèi)流線的對比

圖5 一個周期內(nèi)的阻力系數(shù)對比

3 算例2——低雷諾數(shù)圓柱繞流

3.1 定常計算

首先，以Re=41的定常圓柱繞流為例驗證所采取的預處理方法.計算域為20倍圓柱直徑，網(wǎng)格為91×81的簡單O型網(wǎng)格，流動的馬赫數(shù)為0.001.從圖6可以看到不使用預處理時，耗散過強，流動沒有分離;而采用預處理后，減少了對流格式的數(shù)值耗散，計算的流線和實驗吻合很好，分離區(qū)的大小約為2倍圓柱直徑.

圖6 Re=41的低速圓柱繞流的流線圖

3.2 非定常計算

當圓柱繞流的雷諾數(shù)增加到Re=200時，圓柱后就會出現(xiàn)周期性的渦的脫落，出現(xiàn)“卡門渦街”.首先采用非定常計算，預處理馬赫數(shù)取為0.01，計算網(wǎng)格為156×156，其計算的阻力系數(shù)Cd和升力系數(shù)Cl隨時間的演化過程如圖7所示.計算得到的阻力系數(shù) Cd=1.35±0.039，升力系數(shù)Cl= ±0.64，渦脫落頻率(Strouhal數(shù))St=0.19.本文的計算結(jié)果和其他研究者的計算或?qū)嶒灲Y(jié)果吻合很好，如表1所示.

圖7 阻力、升力系數(shù)隨時間的變化

表1 Re=200的圓柱繞流的非定常計算結(jié)果對比

在采用諧波平衡法計算時，需要事先知道渦的脫落頻率，或者采用參考文獻［13］提出的一種基于梯度類的尋優(yōu)算法來動態(tài)的決定渦的脫落頻率.這里直接利用時間推進法計算得到的渦的脫落頻率作為諧波平衡法中渦脫落的基準頻率.

采用不同階數(shù)諧波計算的阻力、升力系數(shù)如表2所示.隨著所采取的諧波階數(shù)的增加，阻力系數(shù)的均值保持在1.35左右，和時間推進計算結(jié)果基本一致;升力系數(shù)為±0.67，比時間推進計算結(jié)果略大，但在表1其他研究者計算的范圍之內(nèi).

表2 不同階諧波計算得到的阻力、升力系數(shù)

在諧波平衡法計算中，當采用較少階數(shù)諧波時，計算誤差大;采用更多諧波計算時，計算精度雖然會得到改善，但增加了計算時間和內(nèi)存需求.為了分析實際需要的合理的諧波階數(shù)，在圓柱后沿軸向選擇4個監(jiān)控點(A，B，C，D).利用時間推進法計算得到的一個周期內(nèi)守恒變量ρu的變化規(guī)律如圖8所示.可以看到在圓柱后渦交替脫落的A點，其在一個周期內(nèi)的變化劇烈，需要采用較多階數(shù)的諧波來分辨;而隨著渦向下游的不斷傳遞，由于粘性耗散，其ρu的變化趨于平緩，只需要較少階數(shù)的諧波.這也可以從圖9，這4個點ρu的傅里葉分析的頻譜圖中看出，在A點，其前6階諧波的幅值和一階諧波幅值相比，不可忽略，故需要大約6階諧波才能較好的分辨出A點的非定常信號，而在B，C，D點需要的諧波階數(shù)越來越少，在D點，只需3階諧波即可.

圖8 時間推進法計算的一個周期內(nèi)監(jiān)控點ρu的分布

圖9 時間推進法計算的監(jiān)控點的ρu的頻譜

圖10為采用不同階數(shù)諧波計算還原得到的t=0.5T時刻的流場和時間推進法計算流場的對比.采用一階諧波，能夠分辨出渦脫落的基本過程，但是分辨率很差.采用3階諧波對C，D點附近的流場細節(jié)有很好的分辨能力，但對A，B點附近的流場分辨不夠，這和圖9的傅里葉分析結(jié)果匹配.當采用6階諧波計算時，就可以很好的還原出整個渦脫落區(qū)的流場細節(jié).

圖10 t=0.5T時刻渦量的對比

上面的分析表明如果能夠采用自適應的諧波平衡法，即是在流場變化大的區(qū)域采用較多階數(shù)的諧波，在流場變化平緩的區(qū)域采用較少階數(shù)的諧波，這樣與全場采用同樣多高階諧波的計算相比，可以減少計算量和內(nèi)存需求，這也是作者后續(xù)研究的方向之一.

4 結(jié)論

1)采用時間預處理方法，可以用基于可壓縮流方程的計算格式求解低速不可壓流，從而將高、低速流的計算在同一個方程中統(tǒng)一起來.

2)在低速周期性非定常流計算中，帶預處理的諧波平衡方程可以采取較少階數(shù)的諧波計算，還原得到一個周期內(nèi)幾乎準確的非定常流場.表明在計算規(guī)模大，采用時間推進法求解需要漫長的中間過渡態(tài)的周期性非定常流計算中，采用諧波平衡法，可以大大減少計算時間.

3)有必要進一步研究自適應的諧波平衡法，減少諧波平衡法的計算時間和內(nèi)存需求.

References)

［1］ He L.Fourier methods for turbo machinery applications［J］.Progress in Aerospace Science，2010，46(8):329 －341

［2］ Hall K C，Thomas JP，Clark W S.Computation of unsteady nonlinear flows in cascadesusing a harmonic balance technique［J］.AIAA Journal，2002，40(5):879 －886

［3］ Lucia D J，Beran P S，Silva W A.Reduced-order modeling:new approaches for computational physics［J］.Progress in Aerospace Sciences，2004，40(1 －2):51 －117

［4］ Turkel E.Preconditioning techniques in computational fluid dynamics［J］.Annual Review of Fluid Mechanics，1999，31:385 －416

［5］ Briley W R，Taylor L K，Whitfield D L.High-resolution viscous flow simulations at arbitrary Mach number［J］.Journal of Computational Physics，2003，184:79 －105

［6］ Sicot F，Puigt G，Montagnac M.Block-Jacobi implicit algorithms for the time spectralmethod［J］.AIAA Journal，2008，46(12):3080－3089

［7］ Edwards JR.Towards unified CFD simulationsof real fluid flows［R］.AIAA-2001-2524，2001

［8］曹寧，吳頌平.低速流預處理Roe格式中的數(shù)值粘性［J］.北京航空航天大學學報，2010，36(8):904－908 Cao Ning，Wu Songping.Numerical dissipation of Roe's scheme with preconditioning for low-speed flows［J］.Journal of Beijing University of Aeronautics and Astronautics，2010，36(8):904 －908(in Chinese)

［9］ Ghia U，Ghia K N，Shin C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and amultigridmethod［J］.Journal of Computational Physics，1982，48:387 － 411

［10］ Dyke Van，Milton D.An album of fluid motion［M］.Stanford，CA:Parabolic Press，1982

［11］ Chin Hoe Tai，Zhao Yong.Parallel unsteady incompressible viscous flow computations using an unstructured multigrid method［J］.Journal of Computational Physics，2003，192:277 － 311

［12］ Rosenfeld M，Kwak D，Vinokur M.A solution method for the unsteady and incompressible Navier-Stokes equations in generalized coordinate systems［R］.AIAA-88-0718，1988

［13］ McMullen M，Jameson A，Alonso JJ.Application of an on-linear frequency domain solver to the Euler and Navier Stokes equations［R］.AIAA 2002-0120，2002

(編 輯:張 嶸)

Application of harmonic balance method in simulations of low speed unsteady flows

Du Pengcheng Ning Fangfei
(School of Jet Propulsion，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

The harmonic balance method is an effective computational method in simulating time periodic unsteady flows.When using the harmonic balance method based on compressible flow equations to solve low Mach number flows，both the accuracy and convergence of the solution would be deteriorated due to the large dissipation of the convective scheme which is specially designed for compressible flows.To solve this problem，the low speed preconditioning was adopted;therefore，the harmonic balance method based on compressible flow equations could be used to compute low speed periodic unsteady flow directly.Both the time marching calculation and the harmonic balance calculation were performed in simulating the incompressible lid-driven flow and the low Reynolds number vortex shedding cylinder flow.The results show the capability of using the preconditioned harmonic balance equation in simulating low speed periodic unsteady flows，and the unsteady flowfield could be well reconstructed by using only a few harmonics retained in the harmonic balance method.

periodic unsteady flow;harmonic balance method;low speed preconditioning;lid-driven flow;cylinder flow

V 211.3

A

1001-5965(2012)06-0766-06

2011-03-16;網(wǎng)絡出版時間:2012-06-15 15:44

www.cnki.net/kcms/detail/11.2625.V.20120615.1544.031.htm l

國家自然科學基金資助項目(50506001)

杜鵬程(1986－)，男，四川宣漢人，博士生，dupengcheng22@163.com.