基于方、矩形環(huán)管直接數(shù)值模擬結(jié)果的各向異性RANS湍流模型構(gòu)建

郭永濤，徐弘一，2

(1.復旦大學 航空航天系，上海 200433； 2.大連理工大學 工程力學系國際計算力學研究中心 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116023)

歷經(jīng)一個半多世紀的研究，湍流問題仍然是一個懸而未決的經(jīng)典物理學問題.非線性Navier-Stokes(N-S)方程在經(jīng)過雷諾平均運算后，產(chǎn)生了使方程無法封閉的雷諾應力項.由于受到計算機發(fā)展水平的限制，當今工業(yè)界的研究焦點仍在探索如何準確地模型化這些雷諾應力項，而非求助于直接數(shù)值模擬(Direct Numerical Simulation，DNS)技術(shù).主要原因是DNS目前還只能應用于較為簡單的流動幾何構(gòu)型，如槽道流[1]，圓管流[2]及同心圓管流[3]，方、矩形環(huán)管[4]等.而對于復雜構(gòu)型的湍流流動，DNS目前尚無能為力，近期僅能作為一些湍流基礎(chǔ)研究的手段，如湍流的轉(zhuǎn)捩現(xiàn)象及內(nèi)在機理[5-7]、湍流減阻[8-10]等.因此在今后的相當長時間內(nèi)，雷諾平均方法(Reynolds-Averaged Navier-Stokes, RANS)還將是實際工程應用的主要選擇.

1877年，Boussinesq[11]首先根據(jù)湍流與分子運動的類比，提出了雷諾應力和平均應變率在平板流動中的關(guān)聯(lián)關(guān)系，而后Kolmogorov[12](1942年)提出了3維湍流的一般湍流渦粘假設(shè)和本構(gòu)關(guān)系，表達式如下：

(1)

徐弘一[4,13,19-20]模擬所得到的方、矩形環(huán)管流動展現(xiàn)了較為豐富的湍流流動信息，并為本文研究提供了方、矩形環(huán)管DNS數(shù)據(jù)庫.本文進一步根據(jù)徐弘一[4]所提出的各向異性湍流本構(gòu)關(guān)系建立各向異性湍流模型，并利用DNS數(shù)據(jù)庫進行校正和驗證.徐弘一已在文獻中將得到的計算數(shù)據(jù)與其他研究人員的實驗數(shù)據(jù)和DNS數(shù)據(jù)進行對比，驗證了數(shù)據(jù)的可靠性，這也為本文的計算和驗證奠定了基礎(chǔ).

在流體分子運動各向同性的基礎(chǔ)上，Navier和Stokes分別推導出流體力學最為重要的基本方程之一——N-S方程.徐弘一[4，13]的研究表明一般各向異性粘性系數(shù)4階張量在各向同性假設(shè)下，可退化為含有2個獨立標量粘性系數(shù)的各向同性4階張量.由此得到的線性本構(gòu)關(guān)系，即為傳統(tǒng)的牛頓流體本構(gòu)關(guān)系.在仿照該本構(gòu)關(guān)系建立湍流本構(gòu)關(guān)系時，不妨將各向異性湍流應力與各向異性湍流應變率進行關(guān)聯(lián)，而相應的渦粘系數(shù)則應該表述為4階各向異性張量.徐弘一[4，13,19-20]已利用直接數(shù)值模擬的數(shù)據(jù)從應力場分布指出湍動能的各向異性，并利用能譜分析與Kolmogorov的-5/3律比對發(fā)現(xiàn)只有管流中心區(qū)附近的小尺度運動符合各向同性假設(shè).文獻[21-23]運用張量識別定理亦得出張量形式的渦粘系數(shù)應為4階張量.綜上所述，可得出更加一般的適合各向異性湍流運動的本構(gòu)關(guān)系：

(2)

式中：λij為各向異性湍動能；υijkl為張量渦粘系數(shù)；下標i,j,k,l=1,2,3分別表示x,y,z方向；重復下標表示愛因斯坦求和.

由于指標i和j對稱，k和l對稱以及不可壓性條件，渦粘系數(shù)張量獨立分量可以減少到30個[21-23].雖然工作量已大大減少，但其分量仍太多，無法在計算中應用.雖然文獻[21-23]使用廣義逆矩陣求解張量渦粘系數(shù)，但得到的僅為最小模解，很難說是真實的渦粘系數(shù).Carati[24]提出了用2個4階張量來關(guān)聯(lián)雷諾應力和平均速度，但模型過于復雜以至于難以求解.Nieckele[25]近期采取一種非線性渦粘模式來模型化各向異性湍流.本文將分析對稱性較高的方、矩形環(huán)管的充分發(fā)展湍流數(shù)據(jù)，通過幾何對稱性以期減少未知量個數(shù)，并提取主要影響分量進行數(shù)值計算，與DNS數(shù)據(jù)對比以驗證上述采用張量形式渦粘模型表述各向異性湍流運動的猜想.

1 充分發(fā)展湍流方程

假設(shè)在空間某一位置流體的物理量瞬時值等于平均和脈動值之和，并結(jié)合N-S方程可以得到雷諾平均方程：

(3)

(4)

式中：t為時間；p為壓強；Re為雷諾數(shù).

(5)

(6)

(7)

(8)

2 基于直接數(shù)值模擬結(jié)果的渦粘系數(shù)構(gòu)造

2.1 雷諾應力及應變率分布

圖1和圖2展示了直接數(shù)值模擬得到充分發(fā)展湍流的方、矩形環(huán)管的雷諾應力分布.其中基于特征長度h和平均摩擦速度uτ的雷諾數(shù)為Reτ=uτh/ν=200.

圖1 方環(huán)管雷諾應力分布Fig.1 Reynolds stress distribution of square annular duct

圖2 矩形環(huán)管雷諾應力分布Fig.2 Reynolds stress distribution of rectangular annular duct

觀察方、矩形環(huán)管的雷諾應力分布圖，可以看出兩者分布相似，只是方環(huán)管存在旋轉(zhuǎn)90°的對稱性，而矩形環(huán)管僅具有軸對稱特性.

雷諾應力的這種對壁面和角域的選擇特性體現(xiàn)了湍流脈動受宏觀幾何形狀影響較大，也在這些區(qū)域反映出湍流的各向異性這一基本特征.

圖3 平均應變率分布Fig.3 Average strain rate distribution

2.2 渦粘系數(shù)模型

(9)

(10)

(11)

(12)

在用直接數(shù)值模擬數(shù)據(jù)和方程(11)，(12)中的渦粘系數(shù)來構(gòu)造上述湍流應力本構(gòu)關(guān)系時，發(fā)現(xiàn)渦粘系數(shù)全場絕大多數(shù)為正值(見圖4)，僅在管流中心零應變率區(qū)由于應變率在此處接近于0，造成計算結(jié)果出現(xiàn)諸多壞點.在壞點剔除后，所得到的渦粘系數(shù)的恒正性滿足了熱力學第二定律要求.

而展向雷諾切應力與應變率應滿足：

(13)

2.2.1 對稱渦粘系數(shù)模型

對于方形環(huán)管，其幾何上存在著90°旋轉(zhuǎn)對稱性，因此在這里可以假設(shè)張量渦粘系數(shù)也存在這樣的對稱性，即2υ1212=2υ1313=α1和2υ1213=2υ1312=α2.于是可以簡化原方程至只含2個相互獨立的未知渦粘系數(shù)方程：

(14)

(15)

求解以上關(guān)于α1和α2的二元一次方程得到：

(16)

α1和α2作為所關(guān)聯(lián)出的渦粘系數(shù)，因其隱含了對稱假設(shè)，稱之為對稱渦粘模型.

2.2.2 雙渦粘系數(shù)模型

從前面的分析及圖2，圖3可以看出，υ1213與υ1312對雷諾切應力和應變率的關(guān)聯(lián)度不高.在對稱渦粘系數(shù)模型中，可以計算出2個渦粘系數(shù)，見圖4，發(fā)現(xiàn)υ1212遠大于υ1213，因此計算時可以舍去υ1213.若簡化時認為可舍去υ1213和υ1312，則可以分別求得2個渦粘系數(shù)υ1212和υ1313，而且對于一般的矩形管也適用，即2υ1212=β1和2υ1313=β2.同樣得到2個渦粘系數(shù)：

(17)

β1和β2作為舍去了υ1213和υ1312所得到的渦粘系數(shù)，稱之為雙渦粘系數(shù)模型.

2.3 渦粘系數(shù)計算結(jié)果及分析

2.3.1 方環(huán)管渦粘系數(shù)

根據(jù)雙渦粘系數(shù)模型求出的β1和β2分布如圖5所示.

圖4 方環(huán)管對稱渦粘系數(shù)分布云圖Fig.4 Contours of symmetric eddy viscosity distribution of square annular duct

圖5 方環(huán)管雙渦粘系數(shù)分布云圖Fig.5 Contours of double eddy viscosity distribution of square annular duct

2.3.2 矩形環(huán)管渦粘系數(shù)

矩形環(huán)管僅具有軸對稱特性，而不具有方環(huán)管的旋轉(zhuǎn)90°對稱性.分析文獻[20]中圖3(b)和圖4(b)，可發(fā)現(xiàn)主流方向上矩形管在凹角和凸角附近依然可保持一定的對稱性，而凹凸角正是湍流模擬的重要區(qū)域.因此，在利用對稱渦粘系數(shù)模型求解時，不妨強加給矩形環(huán)管90°旋轉(zhuǎn)對稱性進行求解.根據(jù)對稱渦粘系數(shù)模型求出的α1和α2分布如圖6所示.

對矩形環(huán)管采用雙渦粘系數(shù)模型求出的β1和β2分布如圖7所示.

圖6 矩形環(huán)管對稱渦粘系數(shù)分布云圖Fig.6 Contours of symmetric eddy viscosity distribution of rectangular annular duct

圖7 矩形環(huán)管雙渦粘系數(shù)分布云圖Fig.7 Contours of double eddy viscosity distribution of rectangular annular duct

2.3.3 異常數(shù)據(jù)處理及兩種渦粘模型比較

圖8 z=2.955(0≤y≤1)處渦粘系數(shù)分布圖Fig.8 Eddy viscosity distribution at z=2.955(0≤y≤1)

在圖8中，中間3個異常點(2個較大，1個小于0)，明顯與周圍的渦粘系數(shù)成不同的變化趨勢，因此判斷為不合理的渦粘系數(shù)值.運用鄰近的2個合理值進行線性插值以期修正或改善異常的渦粘系數(shù)的影響.

經(jīng)過修正的渦粘系數(shù)分布圖如圖9～圖12，其中方、矩形環(huán)管對稱渦粘系數(shù)的α2修正后數(shù)值相對于α1過小而予以舍去.

在修正的3個異常點中，負渦粘系數(shù)對計算影響較大，可導致計算發(fā)散.這是由于正渦粘系數(shù)表征為湍流引起的能量耗散和動量的傳遞方向——從大渦傳遞至小渦，而負渦粘系數(shù)則意味著小渦將動量輸運至大渦，在簡單剪切流中顯然不合理.2個較大的渦粘系數(shù)經(jīng)過插值近似后使得渦粘系數(shù)分布變得光滑，但由于管流中心區(qū)應變率近似于零，其作用則相對較小.

對于雙渦粘模型得到的系數(shù)，如圖10和圖12，雖經(jīng)過修正，仍然可以看出的β1上下兩邊和β2左右兩邊仍有少量間斷和非光滑過渡.通過觀察應變率的分布情況，y方向切應變對上下兩邊而言數(shù)值較小，如圖3(a);z方向切應變在上下兩邊所轄區(qū)域數(shù)值較小，如圖3(b).故可知渦粘系數(shù)存在間斷區(qū)域?qū)^小的切應變率分布，因而對計算結(jié)果影響較小.

根據(jù)前面的分析，可知間斷出現(xiàn)的區(qū)域?qū)^小的應變率分布，不準確的渦粘系數(shù)因與接近0值的應變率相乘，使計算得到的雷諾應力誤差會相對比較小.因此，光滑后的渦粘系數(shù)分布雖存在若干不光滑點，但對計算結(jié)果影響不大.

圖9 方環(huán)管修正對稱渦粘系數(shù)α1分布云圖Fig.9 Contours of modified symmetric eddy viscosity α1 distribution of square annular duct

圖10 方環(huán)管修正雙渦粘系數(shù)分布云圖Fig.10 Contours of modified double eddy viscosity distribution of square annular duct

圖11 矩形環(huán)管修正對稱渦粘系數(shù)α1分布云圖Fig.11 Contours of modified symmetric eddy viscosity α1 distribution of rectangular annular duct

圖12 矩形環(huán)管修正雙渦粘系數(shù)分布云圖Fig.12 Contours of modified double eddy viscosity distribution of rectangular annular duct

3 雷諾平均方程求解算法及結(jié)果比較

3.1 雷諾平均方程求解算法

由于式(6)，(7)中的雷諾應力項難以用渦粘模型構(gòu)造，在前面的敘述中也從數(shù)值求解的難度上指出了展向面內(nèi)構(gòu)造渦粘系數(shù)的不可靠性.湯泓灝[26]從二次流角度指出右端的雷諾應力項不能用渦粘系數(shù)去構(gòu)造，渦粘模型主要體現(xiàn)了湍流的耗散性，而無法體現(xiàn)湍流的產(chǎn)生，即二次流.Speziale[27-28]指出，非線性的湍流模型可以模擬出二次流現(xiàn)象，但得到的結(jié)果仍不能反映湍流的全部機理[29].

將模型化的湍流本構(gòu)關(guān)系帶入式(6)可得：

(18)

然而對式(18)的直接離散會造成計算精度的損失，不妨將雷諾應力項與擴散項合并以進行整體離散，得到：

(19)

對式(19)進行離散時，對流項采用1階Adams-Bashforth顯式格式，擴散項采用2階Adams-Moulton隱式格式.同樣，渦粘系數(shù)項需和擴散項一起整體離散，離散時需注意到由于υ1212為空間的函數(shù)，不能當作常數(shù)處理.進行整體離散的好處在于可以避免對υ1212進行多次差分運算，從而提高計算精度和穩(wěn)定性.

對于方、矩形環(huán)管，采用雷諾應力作為源項進行計算時，高度還原了直接數(shù)值模擬所得到的充分發(fā)展湍流中的平均速度場數(shù)據(jù).本文所得到的方、矩形環(huán)管展向速度分布如圖13所示.

圖13 環(huán)管展向合速度分布及矢量圖Fig.13 Resultant velocity distribution of annular duct cross-section

圖14 SST k-ω模型得到的方環(huán)管展向合速度 分布及矢量圖Fig.14 Resultant velocity distribution of square annular duct cross-section from SST k-ω model

對于方環(huán)管，用對稱渦粘系數(shù)模型計算所得全場平均相對誤差為2.39%；雙渦粘系數(shù)模型計算所得全場平均相對誤差為1.62%.其中后者相對誤差較小的原因可能在于對流場信息還原比較全面，而前者的缺陷則主要由計算中舍棄了1個渦粘系數(shù)所致.從速度分布對比圖15及圖16可以看出，管流中心區(qū)的絕對誤差較大，這是管流中心區(qū)網(wǎng)格較粗且處于零應變率區(qū)造成的.

圖15 各模型得到的方環(huán)管流向速度分布對比圖 distribution of square annular duct from different model

圖16 方環(huán)管截面不同位置處的速度分布Fig.16 Velocity distribution at different position of square annular duct cross-section

圖17 矩形環(huán)管截面不同位置處的速度分布Fig.17 Velocity distribution at different position of rectangular annular duct cross-section

對于矩形環(huán)管，用對稱渦粘系數(shù)模型計算所得全場平均相對誤差為6.83%；雙渦粘系數(shù)模型計算所得全場平均相對誤差為5.57%.速度分布與方環(huán)管類似.但其誤差要比方環(huán)管要大，且渦粘模型計算結(jié)果均大于DNS結(jié)果，原因在于對于渦粘系數(shù)處理時插值的渦粘系數(shù)精度降低，導致耗散性不夠，但誤差也一般在工程允許范圍之內(nèi).

無論運用于方環(huán)管還是矩形環(huán)管，對稱渦粘系數(shù)模型都比雙渦粘系數(shù)模型計算所得全場平均相對誤差要大.因此雙渦粘系數(shù)模型的計算精度較高，而對稱渦粘系數(shù)模型由于舍去了一個渦粘系數(shù)，計算量會稍低一些.而用SSTk-ω基準模型計算的誤差普遍較大，尤其在角域，這是因為SSTk-ω模型無法在角域預測出二次流，使得速度剖面誤差較大.對比可知，本文發(fā)展的張量渦粘模型要明顯優(yōu)于SSTk-ω模型.

4 結(jié) 語

本文指出了傳統(tǒng)標量渦粘模型的缺陷，雷諾正應力和展向切應力無法用渦粘系數(shù)與應變率進行關(guān)聯(lián).在張量渦粘系數(shù)的假設(shè)下，本文建立了1種簡單的各向異性RANS模型來模型化主流方向的雷諾切應力，并利用方、矩形環(huán)管的直接數(shù)值模擬結(jié)果驗證了模型的適用性，為后續(xù)基于張量渦粘系數(shù)來模型化湍流雷諾應力提供了1個可行方向.計算證明利用湍流模型可大大減少計算時間，完成整個RANS的計算僅需2h的個人電腦CPU時間(Intel(R) Core(TM) i7-4790 CPU@3.60GHz 8核，內(nèi)存32GB)，而相關(guān)DNS計算則需數(shù)月的大型服務器CPU時間.然而，由于目前可靠的DNS數(shù)據(jù)有限，湍流模型構(gòu)建只能應用于少數(shù)流動幾何構(gòu)型簡單的流場.本文僅成功地模型化了主流方向的流動，后續(xù)工作需要開發(fā)其他湍流模型(非渦粘模型)對展向雷諾應力進行模擬.在這方面，近期Ling[32]和Kutz[33]基于大數(shù)據(jù)的深度神經(jīng)網(wǎng)絡學習方法為湍流模型的構(gòu)建提供了1個誘人前景.因此，在本文對DNS結(jié)果中各雷諾應力項的精確驗證基礎(chǔ)上，下一步工作將以DNS大數(shù)據(jù)結(jié)果作為訓練樣本，采用深度神經(jīng)網(wǎng)絡學習的方法來嘗試構(gòu)建湍流模型，以解決雷諾正應力和展向切應力的模型化難題.

致謝： 本研究部分工作在國家超級計算天津中心天河計算機上完成.同時感謝曹博超老師對本文未來工作的建議.