低階格式在大渦模擬計算中的適用性

劉同新， 馬寶峰

（北京航空航天大學流體力學教育部重點實驗室，北京 100191）

低階格式在大渦模擬計算中的適用性

劉同新， 馬寶峰

（北京航空航天大學流體力學教育部重點實驗室，北京 100191）

采用三維Taylor-Green渦作為研究對象，利用工程中常用的低階數值格式，研究格式本身的數值誤差對大渦模擬計算的影響.結果表明：三種數值格式的數值耗散行為都與亞格子模型行為類似，即在小雷諾數下，流場比較光滑時，耗散很小，當雷諾數增加，流動轉捩為湍流，流場梯度增大，耗散顯著增大.對于MUSCL格式和二階有界中心格式，在高雷諾數下，亞格子尺度模型沒有明顯改善計算結果，但也沒有使計算結果惡化.中心格式相比其它兩種格式，數值耗散最小，但是在高雷諾數湍流情況下，中心格式的數值耗散仍然主導了能量的耗散，再添加亞格子模型，計算結果反而變得稍差.對于工程中的低階格式而言，采用中心格式計算大渦模擬是比較好的選擇，而且在計算不存在穩定性問題時，采用不添加亞格子模型的隱式大渦模擬效果更好.

大渦模擬；數值耗散；亞格子模型；低階格式；

0 引言

湍流的準確計算在航空、航天器設計和研制中具有重要意義，因而能夠準確地模擬湍流一直是學術界和工業界追求的目標.在工程湍流計算中以往大都采用基于雷諾平均的湍流模式（RANS），由于湍流模式是對湍流內大小不同尺度的流動結構統一建模，難以得到普適模型，計算準確度欠佳.同時雷諾平均基于時間平均，抹掉了非定常脈動細節，無法計算轉捩、氣動噪聲等非定常問題.直接數值模擬方法（DNS）盡管計算準確，并可獲得流場的全部信息，但計算資源的要求過高，在可以預見的將來難以應用到工程計算中，目前主要用來對低雷諾數、簡單幾何外形的流動做機理研究.

大渦模擬（LES）是目前最有希望取代雷諾平均來進行湍流工程計算的方法［1］.其在相對較低的計算資源下能夠模擬較高的雷諾數運動，同時能夠獲得非定常流場信息.大渦模擬不對時間求平均，而是在空間上對湍流小尺度渦進行濾波，因而流動在時間方向上大渦脈動信息不丟失.濾波掉小尺度渦后由于切斷了湍流從大渦向小渦的級串，因而能量會在截斷的尺度上堆積，產生非物理解或數值不穩定，因而在計算時需要設計一種機制，將累積的能量耗散掉.用來耗散小渦能量的模型稱為亞格子尺度模型（subgrid-scale model，SGS）.亞格子模型是針對小渦建模，而小渦運動更容易趨向均勻各向同性，因而更容易獲得相對普適的模型.

對于亞格子尺度模型，除了本身的建模誤差外，模型提供的物理耗散與計算格式誤差之間干擾是一個需要解決的重要問題.Ghosal［2］通過理論分析指出即使是四階精度的中心格式，其截斷誤差與亞格子應力也是同一量級的.所以為了準確模擬亞格子應力，大渦模擬一般要求采用高精度低耗散的數值格式進行計算或者采用顯示濾波方法［3］濾除網格分辨率內的高波數段模擬不準的流場.高精度格式和顯示濾波方法雖然在學術研究領域，針對簡單幾何外型的流場取得較好的計算結果，但出于計算效率和魯棒性的原因，一直沒有廣泛應用到工程計算中.工程領域采用的仍然是二階到三階精度的低階格式.基于數值格式誤差對亞格子模型的干擾較大，有研究人員提出不顯式建立亞格子尺度模型，而利用數值格式自身的數值耗散作為隱式模型模擬SGS應力，該大渦模擬方法稱為隱式大渦模擬.隱式大渦模擬類似于超聲速計算中的激波捕捉方法，即都是用數值格式自身的耗散來模擬物理問題.

Boris等［4］首次分析了單調格式的截斷誤差，發現其耗散行為與SGS應力很類似，因而提出可利用格式的截斷誤差隱式作為SGS模型的想法，而不再建立顯式SGS模型.但正如Grinstein等［5］所說的，并不是所有的數值格式不加顯式SGS模型，在粗網格上計算都叫做隱式大渦模擬，只有某些特定的數值格式才能夠得到較好的計算結果.多種數值格式已被驗證適合做隱式大渦模擬，可參照相關綜述文章［5］.Garnier等［6］采用均勻各向同性湍流，研究了一系列迎風格式在大渦模擬中的適用性，發現迎風類格式的數值耗散行為與SGS模型類似但不完全一樣，即使是5階WENO格式的數值耗散已足夠大，不建議另外再加顯式SGS模型.而Hickel［7］等進一步指出，利用已有的格式采用逐個嘗試的方法研究隱式大渦模擬有些盲目，他們采用改型方程分析方法，針對格式的截斷誤差性質，專門建立了一種適合做隱式大渦模擬的新格式，針對模型流動，取得了很好的模擬結果.Kravchenko and Moin［8］對渠道流動進行大渦模擬誤差分析，他們利用譜方法進行分析，發現在低階格式下誤差對于亞格子模型具有較大的影響，但仍然能夠得到較好的計算結果.楊小龍和符松［9］采用pade格式，研究了低階到高階格式數值誤差對大渦模擬的影響，結果表明亞格子模型能夠緩解數值誤差對計算結果的影響，采用低精度格式也可能夠得到較好計算結果.謝志剛等［10］也采用二階格式得到了較為滿意的計算結果.

綜上，目前對于數值耗散對大渦模擬計算結果的影響，結論尚不一致，不同研究者采用不同的計算格式和方法可能會得到不同的結論，因而需要開展更廣泛的研究，對不同的數值格式特別是工程中常用的低階格式在大渦模擬中的適用性，采用不同的方法進一步進行驗證，從而為工程計算提供參考.本文采用一種后驗的方法（a posteriori），以Taylor-Green渦為計算對象，對工程中常用的低階數值格式的數值耗散對大渦模擬的影響進行評估.

1 數值方法

1.1 主控方程及亞格子尺度模型

本研究的所有計算基于FluentTM14.5代碼，該代碼目前在工業界應用較多.

計算基于不可壓N-S方程，大渦模擬需要對方程進行濾波處理，濾波后的方程為

大渦模擬的亞格子尺度建模既是針對ij進行，本研究的亞格子尺度模型選用動態Smagorinsky模型，這是目前大渦模擬計算中的主流模型，在計算效率和準確性方面較適合于工程計算.動態Smagorinsky模型是將經典Smagorinsky模型中的系數Cs在計算過程中根據大渦信息自動計算［11－12］.

經典Smagorinsky模型基于Boussinesq渦粘假設

動態亞格子建模過程是一種方法，而不是單一的模型，該過程可應用于任何常系數亞格子模型，將其變為動態模型.應用于經典Smagorinsky模型則構成動態Smagorinsky模型.該方法的基本思路是對流動方程進行二次濾波，根據第二次濾波及大渦模擬第一次濾波信息自動計算渦粘系數，從而保證亞格子模型可適應更廣泛的流動計算.

1.2 離散格式

數值離散基于非結構網格有限體積法，不可壓流動方程采用分數步方法（fractional step）［13］求解.相比于SIMPLE類方法，分數步方法不將算子分裂誤差通過迭代降到接近零，而是將分裂誤差降到與時間截斷誤差為同一量級.這樣能夠大大減少在每個時間步內推進的迭代次數，從而在非定常計算中能夠提高計算效率.該算法時間截斷誤差和算子分裂誤差都由時間步長確定，因而時間步長不宜太大.本研究時間離散采用隱式二階格式，具有二階精度，時間步長取為0.001.

本計算空間對流項采用三種離散格式，包括二階中心格式、三階MUSCL格式和二階有界中心格式.粘性項離散采用二階中心格式.

Fluent程序中的有限體積法基于格心法，即計算變量位于各網格單元的幾何中心.不同的空間格式體現為，由相鄰單元中心變量值，計算界面上的值所用的差值或重構方法.

二階中心格式：二階中心格式具有較少的數值耗散常被用在大渦模擬中.該格式利用界面兩側單元中心的值計算界面變量值.對于界面流動變量，基于二階中心格式的界面重構公式為

其中0和1代表了相鄰兩個單元，其中中間的公用面為f，?φr，0和?φr，1是位于單元0和1中心的變量的梯度，其中r為單元中心到面中心的向量.單元中心的變量計算采用Taylor展開法［14］.

因為二階中心格式奇偶失聯，計算時會出現非物理振蕩，因而需要做校正（deferred correction），φf，up為低階迎風格式計算的結果，“old”為上一時間步的結果，當計算收斂后，結果收斂到二階中心格式.

三階MUSCL格式：MUSCL格式是通過中心格式和二階迎風格式混合構造的格式，最早由Van Leer提出［15］.MUSCL格式相比二階迎風格式，具有更小的數值耗散，對網格類型適應性更強.

其中φf，up＝φ＋?φ·r，為二階迎風格式.

二階有界中心格式：由于純中心格式數值耗散較小，即使采用校正，在梯度大的地方計算結果中仍會出現非物理振蕩.而有界中心格式的數值耗散大于純中心格式，但比中心格式魯棒性好，一般不會出現非物理振蕩.有界中心差分格式是由變量歸一化方法（NVD）［16］并在對流有界準則（CBC）［17］約束下構造的復合格式，該復合格式由中心格式、中心和二階迎風混合格式、一階迎風格式復合而成.

1.3 Taylor-Green渦

Taylor-Green渦（TGV）是研究湍流轉捩、級串和耗散的一種重要的模型流動［18］.被廣泛用來驗證大渦模擬的計算正確性.相比均勻各向同性湍流，TG渦包含的流動機制更為豐富，計算難度也更大.

TG渦初始為一光滑的自由周期渦，可視為層流，開始演化后，在非線性作用下，光滑渦結構失穩，轉捩為湍流，并不斷級串為更小的渦，同時隨著時間的演化，能量不斷耗散.

本研究選用三維Taylor-Green渦（TGV）進行計算，分析湍流發展過程中的能量耗散.其初始流場為

初始渦量場可通過ω＝?×u得到.

計算域為邊長為2π的正方體，網格點為64×64×64，正方體對應的兩個面形成周期性邊界條件.根據TG渦直接數值模擬的結果，該網格尺度可保證在所計算的雷諾數下，網格對渦系的截斷波數處于慣性區，從而能夠保證大渦模擬的網格尺度是足夠的.取特征速度為1，特征尺度為1，定義雷諾數為1／ν.計算過程中，雷諾數通過分子粘性系數來改變.

TG渦系統是一個自耗散系統，在演化過程中能量不斷耗散.在數值模擬的情況下，能量耗散的來源分為流體分子粘性耗散，亞格子模型耗散和數值耗散.

分子粘性耗散率 εmol＝2ν＜SijSij＞，＜·＞代表體積平均；

亞格子模型耗散率 εSGS＝2νSGS＜SijSij＞，

其中ν，νSGS，Sij分別分子粘性系數，亞格子湍流渦粘系數和應變率張量，νSGS＝μSGS／ρ．分子耗散率和亞格子耗散率一起稱為物理耗散率εphy＝εmol＋εSGS．

2 計算結果及分析

圖1顯示TG渦在不同時刻的演化結構，初始為規則的周期多渦結構，隨時間演化，逐漸變為小尺度旋渦.圖中渦面是利用渦的λ2判別法［19］得到，其灰度代表了動能的大小.

圖2給出了不同時間步長對計算結果的影響，可見所取時間步長已足夠小，改變時間步長后總能量隨時間的衰減曲線完全重合，本文的以下計算時間步長都取為0.001.圖3給出了不同雷諾數下總平均動能隨時間的演化.在初始動能一樣的情況下，雷諾數越大，耗散越慢，湍流動能在慣性作用下會傳遞到更小尺度的渦上去.

通過與DNS結果的對比，可以比較各種差分格式以及亞格子模型的計算效果，本研究選用BRACHETT［20］的TG渦直接數值模擬結果作為對比數據.BRACHETT［20］的關于TG渦的DNS數據是該方向的經典數據，由譜方法計算，最初計算時采用2563網格，約十年后又基于8643網格做了校核［21］，但未發現偏差.長期以來一直作為TG渦模擬的對比數據，即使最新DNS計算也與其復合得很好［22］.

圖4給出了Re＝100情況下，三種數值格式計算的TG渦系統能量耗散率隨時間的演化.給出了物理耗散率和總耗散率在有無亞格子模型情況下的變化規律，并給出的DNS數據作為對比.可見5條曲線基本上重合到一起，計算結果與DNS數據符合很好，無論加不加亞格子模型，流動的物理耗散率都等于總耗散率，即數值耗散很小.實際上由于在該雷諾數下，流動演化可看成層流演化，耗散過程主要是分子粘性耗散，亞格子耗散可以忽略不計.

圖5給出了Re＝400情況下的計算結果.對于每一種計算格式的計算結果而言，此時物理耗散率和總耗散率都出現差異，特別是在耗散率峰值附近（t＝6 s），差異顯著，表明數值耗散此時對計算結果有明顯的影響.對比三種數值格式，MUSCL格式數值耗散的影響最大，t＝4 s以后，能量總耗散率與物理耗散曲線之間明顯分散（圖5（a）），但是加不加SGS模型對總耗散率和物理耗散率影響不大.中心差分格式和有界中心格式數值耗散的影響相對較小.考慮施加亞格子尺度模型和不加模型兩種情況，采用中心格式計算的物理耗散值與總耗散率的差值是最小的，表明中心格式的數值耗散相對較小，與DNS的結果也比較接近.

當Re＝3 000時，數值耗散的影響顯著增加，如圖6所示.對于MUSCL格式和有界中心格式而言，施加亞格子模型與否，總能量耗散率相差不大.對比DNS數據，亞格子模型僅在局部時間段對耗散率有所改善，但并不明顯.而對于中心格式，不加SGS模型時，總能量耗散率與DNS結果反而更接近，盡管此時耗散仍然偏大，耗散率峰值出現的時間提前，但是峰值的大小與DNS結果很接近.上述結果表明對于三種數值格式而言，在該雷諾數下，數值耗散基本主控了耗散率的量值，且大于亞格子尺度模型提供的耗散.基于這一點可以看出，對于所采用三種低階格式而言，利用格式本身的數值耗散來作為亞格子模型，構成隱式大渦模擬是更為實用的技術，而再顯式添加亞格子模型的效果并不大，甚至會使計算結果變差（比如中心格式）.

對于三種數值格式而言，中心格式的總耗散率與DNS結果最為接近，特別是不加SGS模型時.更為重要的，采用中心格式時，在添加亞格子模型的情況下，中心格式的物理耗散占總耗散率的比重最大（圖6（c）），有界中心格式其次，MUSCL格式最小.這表明中心格式的數值耗散對大渦模擬計算的影響最小.不加亞格子模型時，由于計算物理耗散率時少了亞格子耗散率，所以物理耗散率量值很小.

需要注意，無論加不加亞格子模型，以及采用三種之中的任一數值格式，計算得到的耗散率相比DNS結果在大部分區域都偏大，而在峰值附近又偏小.說明在高雷諾數下還沒有做到準確模擬湍流的耗散.根據已有的研究結果［7］，本文用到動態Smagorinsky模型本身的耗散率相比DNS結果本來就偏大，即該亞格子模型本身也存在建模誤差.因而即使采用高精度低耗散數值格式，顯示添加動態Smagorinsky亞格子模型，計算的耗散率與圖6（c）中計算的總耗散率在趨勢上也是類似的，可參見文獻［7］中四階中心格式的計算結果.

3 結論

以Taylor-Green渦的時間演化為計算對象，研究了三階MUSCL格式、二階有界中心格式和二階中心格式的數值誤差對大渦模擬計算結果的影響.得到如下結論：

1）三種數值格式的數值耗散行為都與亞格子模型行為類似，即在小雷諾數下，流場比較光滑時，耗散很小，當雷諾數增加，流動轉捩為湍流，流場梯度增大，耗散顯著增大.與直接數值模擬結果相比較，對于添加和不添加亞格子模型兩種情況能量耗散都偏大.

2）對于MUSCL格式和二階有界中心格式，在高雷諾數下，亞格子尺度模型沒有明顯改善計算結果，但也沒有使計算結果惡化.中心格式相比其它兩種格式，數值耗散最小，但是在高雷諾數湍流情況下，中心格式的數值耗散仍然主導了能量的耗散，再添加亞格子模型，計算結果反而變得稍差.

3）對于工程中的低階格式而言，采用中心格式計算大渦模擬是比較好的選擇，而且在計算不存在穩定性問題的前提下，采用不添加亞格子模型的隱式大渦模擬效果更好.

［1］張兆順，崔桂香，許春曉.湍流大渦數值模擬的理論和應用［M］.北京：清華大學出版社，2008.

［2］Ghosal S.An analysis of numerical errors in large-eddy simulations of turbulence［J］.Journal of Computational Physics，1996，125（1）：187－206.

［3］Gullbrand J，Chow F K.The effect of numerical errors and turbulence models in large-eddy simulations of channel flow，with and without explicit filtering［J］.Journal of Fluid Mechanics，2003，495：323－341.

［4］Boris J P，Grinstein F F，Oran E S，Kolbe R L.New insights into large eddy simulation［J］.Fluid Dynamics Research，1992，10：199－228.

［5］Grinstein F，Margolin L，Rider W.Implicit large eddy simulation：Computing turbulent flow dynamics［M］.Cambridge，UK：Cambridge University Press，2007.

［6］Garnier E，Mossi M，Sagaut P，Comte P，Deville M.On the use of shock capturing schemes for large-eddy simulation［J］. Journal of Computational Physics，1999，153：273－311.

［7］Hickel S，Adams N A，Domaradzki J A.An adaptive local deconvolution method for implicit LES［J］.Journal of Computational Physics，2006，213：413－436.

［8］Kravchenko A G，Moin P.On the effect of numerical errors in large eddy simulations of turbulent flows［J］.Journal of Computational Physics，1997，131（2）：310－322.

［9］楊小龍，符松.直接數值模擬／大渦模擬中數值誤差影響的研究［J］.應用數學和力學，2008，29（7）.

［10］謝志剛，許春曉，崔桂香，王志石.方柱繞流大渦模擬［J］.計算物理，2007，24（2）：171－180.

［11］Germane M，Piomelli U，Moin P，Cabot W H.A dynamic subgrid-scale eddy viscosity model［J］.Physics of Fluids，1991，（4）：633－635.

［12］Lilly D K.A proposed modification of the germano subgrid-scale closure model［J］.Physics of Fluids，1992，（4）：633－635.

［13］Glaz H M，Bell J B，Colella P.An analysis of the fractional-step method［J］.Journal of Computational Physics，1993，108：51－58.

［14］Barth T J，Jespersen D.The design and application of upwind schemes on unstructured meshes［C］.AIAA 27th Aerospace Sciences Meeting，AIAA－89－0366，Reno，Nevada，1989.

［15］Van Leer B.Toward the ultimate conservative difference scheme.IV.A second order sequel to Godunov's method［J］. Journal of Computational Physics，1979，32：101－136.

［16］Leonard B P.The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection［J］.Computer Methods in Applied Mechanics and Engineering，1991，88：17－74.

［17］Gaskell P H，Lau A K C.Curvative-compensated convective transport：SMART，a new boundedness-preserving transport algorithm［J］.International Journal for Numerical Methods in Fluids，1988，8：617－641.

［18］Taylor G I，Green A E.Mechanism of the production of small eddies from large ones［C］∥Proceedings of the Royal Society of London，Series A，Mathematical and Physical Sciences，1937，158：499－521.

［19］Jeong J，Hussain F M.On the identification of a vortex［J］.Journal of Fluid Mechanics，1991，285：69－94.

［20］Brachet M E，Meiron D I，Orszag S A，Nickel B G，Morf R H，Frisch U.Small-scale structure of the Taylor-Green vortex ［J］.Journal of Fluid Mechanics，1983，130：411－452.

［21］Brachet M，Meneguzzi M，Vincent A.Politano H，Sulem P-L.Numerical evidence of smooth self-similar dynamics and possibility of subsequent collapse for three-dimensional ideal flow［J］.Physics of Fluids，1992，（4）：2845－2854.

［22］Chapelier J-B，Plata M D L L.Renac F.Inviscid and viscous simulations of the Taylor-Green vortex flow using a modal discontinuous Galerkin approach［A］.AIAA－2012－3073，2012.

Suitability of Low-order Numerical Schemes in Large Eddy Simulations

LIU Tongxin，MA Baofeng
（Ministry-of-Education Key Laboratory of Fluid Mechanics，Beihang University，Beijing 100191，China）

Several low-order numerical schemes were evaluated for their suitability in large-eddy simulations based on 3D Taylor-Green vortex，with and without a subgrid-scale model.It shows that dissipation characteristics of three numerical schemes used are similar to a subgrid-scale model.At lower Reynolds numbers，flow fields are relatively smooth，numerical dissipation is lower；As Reynolds number increasing，transition to turbulence occurs and numerical dissipation grows greatly.For MUSCL and bounded centered schemes，subgrid-scale model has a little influence on results at high Reynolds numbers.The second-order central scheme exhibits lower dissipation，but at higher Reynolds numbers numerical dissipation still dominates total energy dissipation of flow.With an explicit subgrid-scale model the results even become worse.Therefore，for large eddy simulations in engineering，the second order central scheme is suitable，particularly without an explicit subgrid-scale model.

large eddy simulation；numerical dissipation；subgrid-scale model；low-order numerical schemes

date： 2013－07－01；Revised date： 2013－10－02

O368

A

1001-246X（2014）03-0307-07

2013－07－01；

2013－10－02

國家自然科學基金（11272033）資助項目

劉同新（1987－），男，山東，碩士生，研究方向：非定常流動的大渦模擬研究，E-mail：bf-ma＠buaa.edu.cn