湍流模型源項(xiàng)的數(shù)值處理方法對計(jì)算結(jié)果的影響研究

王翔宇，李棟

（西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710072）

湍流模型源項(xiàng)的數(shù)值處理方法對計(jì)算結(jié)果的影響研究

王翔宇，李棟

（西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710072）

渦粘性湍流模型廣泛應(yīng)用于各種工程湍流問題的計(jì)算。對于大多數(shù)湍流模型，在湍流控制方程的右端會出現(xiàn)源項(xiàng)，其剛性給數(shù)值計(jì)算帶來很大的影響。從源項(xiàng)的物理意義出發(fā)分析了源項(xiàng)導(dǎo)致求解困難的原因，并在通用的求解湍流模型方程的點(diǎn)隱法基礎(chǔ)上，以源項(xiàng)弱剛性的S-A湍流模型和源項(xiàng)強(qiáng)剛性的k-ωSST湍流模型為例，論述了生成項(xiàng)和耗散項(xiàng)之間的平衡關(guān)系對數(shù)值模擬的重要意義，給出了針對不同源項(xiàng)的具體處理方法。對RAE2882翼型跨聲速流動算例的模擬結(jié)果表明，這些處理方法有效的提高了控制方程組的計(jì)算穩(wěn)定性。

NS方程；源項(xiàng)；剛性；S-A湍流模型；k-ωSST湍流模型；數(shù)值模擬

0 引 言

當(dāng)前，基于渦粘假設(shè)的湍流模型在湍流工程計(jì)算中得到了廣泛的應(yīng)用。與一般NS控制方程相比，渦粘湍流模型方程除了時(shí)間導(dǎo)數(shù)項(xiàng)、對流項(xiàng)和擴(kuò)散項(xiàng)可與之相互對應(yīng)外，最顯著的區(qū)別在于方程中增加了源項(xiàng)。源項(xiàng)的存在一方面可以集中反映了湍流在空間上的尺度多重性和時(shí)間上的高頻脈動性［1］，但另一方面也導(dǎo)致方程數(shù)學(xué)性態(tài)的復(fù)雜變化而產(chǎn)生收斂困難、迭代不穩(wěn)定等一系列問題，甚至可以說源項(xiàng)處理的好壞直接影響著最終對湍流模擬的效果［2］。

基于此，本文首先對湍流模型方程源項(xiàng)數(shù)值求解面臨的數(shù)值剛性問題進(jìn)行討論，其次結(jié)合當(dāng)前應(yīng)用最為廣泛的S-A一方程模型和k-ωSST兩方程模型，較為詳細(xì)的分析了源項(xiàng)的影響，比較了不同的處理手段并以RAE2822翼型跨聲速數(shù)值模擬為算例對結(jié)果進(jìn)行驗(yàn)證。

1 源項(xiàng)的數(shù)值特性分析

一般所指的湍流模型方程源項(xiàng)包括生成項(xiàng)和耗散項(xiàng)兩部分，每項(xiàng)都耦合了時(shí)間特性和空間特性變量，模擬了湍流量生成和消失的過程，在這個(gè)過程中平均流動量演化過程和湍流量演化過程聯(lián)系在一起，對流動計(jì)算的結(jié)果起到了很大的影響［3］。

理論分析指出［2］，湍流方程源項(xiàng)使得流動變量在時(shí)間尺度和空間尺度上發(fā)生了劇烈的變化，而這個(gè)尺度往往遠(yuǎn)小于NS方程的離散尺度，尤其是時(shí)間尺度隨著源項(xiàng)非線性程度的加大會進(jìn)一步減小，即顯著增大了控制方程組的剛性。此時(shí)方程很容易喪失正定性，湍流變量無限制增長，源項(xiàng)大范圍為負(fù)，運(yùn)算很快發(fā)散。從方程穩(wěn)定性上考慮，這些就要求在確定時(shí)間步長時(shí)要分別分析對流項(xiàng)、擴(kuò)散項(xiàng)和源項(xiàng)的當(dāng)?shù)夭介L取其最小值，同時(shí)網(wǎng)格的長寬比、稠密度等也要有更為嚴(yán)格的限制。但這些措施無疑會大大降低控制方程組的收斂速度，而且根據(jù)文獻(xiàn)［4］，很多時(shí)候即便如此也會產(chǎn)生嚴(yán)重的收斂震蕩和非物理解。針對這一問題比較通用的做法是將源項(xiàng)做點(diǎn)隱式處理［4-5］。對于某個(gè)單元i，離散后的方程為：

其中Ω表示這個(gè)單元的體積，W表示控制變量，F(xiàn)kΔSk表示某個(gè)面的各個(gè)通量和，Q表示源項(xiàng)。

可以看出，此時(shí)源項(xiàng)為新時(shí)刻（n＋1）的值，根據(jù)n時(shí)刻的結(jié)果有：

代入式（1）后整理，最終公式為：

為提高控制方程穩(wěn)定性，減少求解時(shí)對CEL數(shù)的依賴，也可用LU-SGS方法隱式求解式（3）：

但由于源項(xiàng)的高度非線性，很多時(shí)候直接求解其雅可比矩陣?Q／?W難度很大，為此可對源項(xiàng)中的生成項(xiàng)和耗散項(xiàng)分別處理，有選擇的隱式求解［6］或者僅對源項(xiàng)的負(fù)部求解雅可比矩陣［6-7］。

另外，生成項(xiàng)和耗散項(xiàng)物理意義不同，這也產(chǎn)生了一個(gè)源項(xiàng)自身平衡的問題。對于許多湍流模型方程，生成項(xiàng)與渦量的模成正比，直觀來看剪切越強(qiáng)，湍流越強(qiáng)，生成項(xiàng)越大。也就是說生成項(xiàng)是湍流強(qiáng)弱的“正尺度”標(biāo)量，其大小密切關(guān)系著湍動能量的注入，是湍流模型方程的主導(dǎo)項(xiàng)，因而對其正定性有著有更為嚴(yán)格的限制。這一點(diǎn)也與零方程模型僅靠正定的生成項(xiàng)自身就能較好模擬附著流動相吻合［2］。對于耗散項(xiàng)，湍流粘性系數(shù)越大，耗散越大，對生成項(xiàng)的限制也就越大。為了在邊界層區(qū)域得到較大的耗散，S-A一方程模型和k-ωSST兩方程模型都引入了壁面距離函數(shù)［8］。但由于壁面附近網(wǎng)格的間距往往非常小，很小的數(shù)值計(jì)算誤差可能被放大到超過生成項(xiàng)自身產(chǎn)生了違反物理本質(zhì)的負(fù)源擴(kuò)散（S-A模型），也可能使得生成項(xiàng)不受耗散項(xiàng)制約增大湍動能無限增長（k-ωSST模型）——這些都將導(dǎo)致數(shù)值模擬的不穩(wěn)定?？紤]到這些問題，后文將對S-A一方程模型和SSTk-ω兩方程模型的生成項(xiàng)和耗散項(xiàng)平衡問題進(jìn)行具體討論并給出相應(yīng)處理辦法。

2 S-A湍流模型的源項(xiàng)處理與算例分析

S-A湍流模型直接假設(shè)其滿足流場中的帶源項(xiàng)的對流-擴(kuò)散標(biāo)量方程，其中包含了大量的經(jīng)驗(yàn)常數(shù)、經(jīng)驗(yàn)函數(shù)。該模型對航空領(lǐng)域的復(fù)雜流動模擬取得了很好的效果。全湍流S-A模型方程可以寫為：

方程中常數(shù)的具體取值可參考文獻(xiàn)［9］。

盡管S-A源項(xiàng)產(chǎn)生的附加剛性較小，很多時(shí)候仍然出現(xiàn)了收斂困難的情況。進(jìn)一步分析，一方面由于耗散項(xiàng)的二次函數(shù)保證其始終為負(fù)，不會出現(xiàn)在迭代計(jì)算中生成項(xiàng)缺乏限制過大產(chǎn)生的湍流量無限增長的情況。但另一方面，壁面處工作變量＝0是SA模型方程在實(shí)際應(yīng)用中給定的邊界條件，當(dāng)壁面網(wǎng)格非常密時(shí)邊界層單元的耗散量較大，很小的（可能非常接近0）無法保證其線性的生成項(xiàng)大于耗散項(xiàng)，導(dǎo)致出現(xiàn)負(fù)源項(xiàng)——這是正常的。但由于原始的S-A模型方程缺乏負(fù)源調(diào)節(jié)機(jī)制，即使只有一個(gè)單元出現(xiàn)負(fù)源，周圍單元正的湍動能量會向其擴(kuò)散導(dǎo)致自身凈生成的減少引起負(fù)源擴(kuò)散，最終使得全場運(yùn)算的發(fā)散。一個(gè)簡單修正的辦法就是在計(jì)算中一旦出現(xiàn)負(fù)的時(shí)（負(fù)源已經(jīng)產(chǎn)生并擴(kuò)散），強(qiáng)制賦予一個(gè)很小的正數(shù)，相當(dāng)于在后續(xù)的迭代中人為注入湍動能量，防止負(fù)源的擴(kuò)散。

根據(jù)上面的討論，下面對經(jīng)典的RAE2882翼型跨聲速流動算例進(jìn)行計(jì)算。首先采用CEL3D網(wǎng)站發(fā)布的RAE2822標(biāo)準(zhǔn)網(wǎng)格，網(wǎng)格點(diǎn)分布為338×65，遠(yuǎn)場約為20倍弦長，第一層網(wǎng)格點(diǎn)到壁面距離約為1×10－5數(shù)量級。計(jì)算條件為［10］：

空間離散為三階MUSCL的Vanleer格式，數(shù)值方法分為以下四種：

solver1：5步Runge-kutta方法（CEL＝0．3）＋源項(xiàng)點(diǎn)隱式計(jì)算；

solver2：LU-SGS方法（CEL＝5．0）＋源項(xiàng)顯式計(jì)算；

solver3：LU-SGS方法（CEL＝5．0）＋源項(xiàng)點(diǎn)隱式計(jì)算；

solver4：LU-SGS方法（CEL＝5．0）＋源項(xiàng)點(diǎn)隱式計(jì)算＋修正。

圖1 RAE2822翼型計(jì)算標(biāo)準(zhǔn)網(wǎng)格Fig.1 Standard grid of RAE2822

圖2 殘值隨迭代變化曲線（S-A模型，標(biāo)準(zhǔn)網(wǎng)格）Fig.2 Residual VS number of iterations（S-A，standard grid）

圖3 翼型表面壓力系數(shù)分布曲線（S-A模型，標(biāo)準(zhǔn)網(wǎng)格）Fig.3 Distribution of pressure coefficient（S-A，standard grid）

從計(jì)算結(jié)果看出，CEL數(shù)取0.3時(shí)顯式Rungekutta方法仍不足以克服近似點(diǎn)隱式處理之后源項(xiàng)的殘余剛性，計(jì)算很快發(fā)散，若一定采用這種方法只能進(jìn)一步降低CEL數(shù)。由于LU-SGS方法自身的穩(wěn)定性，即使源項(xiàng)顯式計(jì)算（迭代方程中略去?Q／?W項(xiàng)，源項(xiàng)顯式計(jì)算在殘差中）也能保證不發(fā)散，但殘差高位震蕩無法收斂，結(jié)果偏差很大。而solver3和 solver4的計(jì)算結(jié)果完全一樣，壓力分布也與實(shí)驗(yàn)數(shù)據(jù)相差很小。這一方面說明了這種“雙隱式”處理很好的模擬了流場，另一方面也說明在當(dāng)前計(jì)算條件下修正并未發(fā)揮作用，即并未出現(xiàn)負(fù)源擴(kuò)散的現(xiàn)象。

考慮到很多分離問題要求對流動變量細(xì)微精確的捕捉，這時(shí)標(biāo)準(zhǔn)網(wǎng)格點(diǎn)分布顯得較為稀疏。接下來對網(wǎng)格加密重新計(jì)算。新網(wǎng)格分布點(diǎn)為368×97，第一層網(wǎng)格點(diǎn)到壁面距離減小到1×10－6數(shù)量級，遠(yuǎn)場仍保持20倍弦長，其余計(jì)算條件都不發(fā)生變化，計(jì)算結(jié)果顯示除slover4外其余三組均無法收斂且其壓力分布與標(biāo)準(zhǔn)網(wǎng)格幾乎一致。進(jìn)一步，圖7給出了距離翼型壁面最近的3層網(wǎng)格單元平均源項(xiàng)值的隨迭代變化的情況。起初solver3和solver4每步平均源項(xiàng)值相同——其中可能有個(gè)別負(fù)源，但總體維持一個(gè)穩(wěn)定狀態(tài)。394步開始可能由于誤差的擾動使得solver3算例中壁面附近大量產(chǎn)生負(fù)源并急劇擴(kuò)散到全場，如圖8所示，之后10步即導(dǎo)致運(yùn)算發(fā)散。而加入修正的solver4算例不存在類似問題，驗(yàn)證了修正可以很好處理S-A模型在邊界層密網(wǎng)格可能產(chǎn)生較大的源項(xiàng)剛性。同時(shí)這個(gè)算例也說明一味加密網(wǎng)格并不能簡單促進(jìn)控制方程對劇烈變化的湍動變量的捕捉，必須綜合考慮各個(gè)因素對控制方程組的影響。

圖4 加密后的RAE2822網(wǎng)格Fig.4 Fine grid of RAE2822

圖5 殘值隨迭代變化曲線（S-A模型，加密網(wǎng)格）Fig.5 Residual VS number of iterations（S-A，fine grid）

圖6 翼型表面壓力系數(shù)分布曲線（S-A模型，加密網(wǎng)格）Fig.6 Distribution ofpressure coefficient（S-A，fine grid）

圖7 翼型附近平均源項(xiàng)值變化曲線（S-A模型）Fig.7 Average of source term near the airfoil VS number of iterations（S-A）

圖8 solver3負(fù)源單元占總單元數(shù)百分比變化曲線Fig.8 Ratio of negative source term to the total vs. number of iterations in solver3（S-A）

3 k-ωSST湍流模型的源項(xiàng)處理與算例分析

k-ωSST兩方程模型是邊界層內(nèi)k-ω模型和遠(yuǎn)場k-ε模型的混合模型。該模型兼顧了兩方面的優(yōu)點(diǎn)：k-ω模型中的壁面衰減函數(shù)，對壁面附近湍流粘性系數(shù)預(yù)測較好；k-ε模型對來流條件不敏感，預(yù)測遠(yuǎn)場流動更為合適。其模型方程可以寫為：

方程中常數(shù)的具體取值可參考文獻(xiàn)［11］。

相比S-A模型方程，k-ωSST模型方程源項(xiàng)要復(fù)雜得多。生成項(xiàng)和耗散項(xiàng)均為湍動的高度非線性函數(shù)，在ω方程中還增加了包含兩個(gè)湍動變量導(dǎo)數(shù)的交叉輸運(yùn)項(xiàng)，這些不但大大提高了控制方程組的剛性，還導(dǎo)致精確求解源項(xiàng)整體的雅可比矩陣變得極為困難。為此需要將源項(xiàng)各部分分開考慮有針對性的進(jìn)行處理。

類似對S-A模型方程的分析，由于k-ωSST模型的固壁邊界條件為

其中d1為距壁面最近一層網(wǎng)格到壁面的距離。一般認(rèn)為［12］不采用壁面函數(shù)的k-ωSST模型壁面最近網(wǎng)格單元當(dāng)?shù)乩字Z數(shù)y＋＜1，壁面網(wǎng)格很密，ω很大，直接導(dǎo)致模型方程耗散項(xiàng)很大，帶來了很大的剛性。因此對于耗散項(xiàng)有必要采用下面的點(diǎn)隱式方法進(jìn)行處理：

在生成項(xiàng)中，

與S-A方程不同，即使壁面附近出現(xiàn)某個(gè)k＜0也不會影響生成項(xiàng)τijSij的正定性，相反卻可能出現(xiàn)湍動能過大的情況。一旦開始迭代時(shí)某單元出現(xiàn)負(fù)的湍動能k，雷諾應(yīng)力τij偏大，生成項(xiàng)偏大，在兩方程模型中邊界層剪應(yīng)力同生成項(xiàng)和耗散項(xiàng)之比的平方根成正比［7］導(dǎo)致剪應(yīng)力偏大，一方面剪應(yīng)力與湍動能的正相關(guān)性會很快將負(fù)的湍動能抹平，另一方面偏大的剪應(yīng)力使得渦量偏大，也就是說生成項(xiàng)會繼續(xù)增大，成為一個(gè)“正反饋”過程。而耗散項(xiàng)卻沒有一個(gè)類似的“機(jī)制”，最終湍動能無限增大計(jì)算崩潰。這也正是前面討論的生成項(xiàng)和耗散項(xiàng)平衡問題，二者的差別不應(yīng)超過一個(gè)數(shù)量級，否則不但影響邊界層流動的模擬能力，更會帶來計(jì)算穩(wěn)定性等嚴(yán)重問題。為此可以將生成項(xiàng)修正為：

即k方程的生成項(xiàng)無論如何不超過其耗散項(xiàng)模的10倍。對于交叉輸運(yùn)項(xiàng)，其自身很大程度上也表示了生成和耗散之間的某種摻混關(guān)系，生成項(xiàng)修正后可直接顯式計(jì)算。

下面采用RAE2822翼型標(biāo)準(zhǔn)網(wǎng)格對SSTk-ω模型進(jìn)行算例驗(yàn)證，計(jì)算條件與S-A模型相同。數(shù)值方法分為下面三種：

solver1：LU-SGS方法（CEL＝5．0）＋耗散項(xiàng)點(diǎn)隱式式計(jì)算；

solver2：LU-SGS方法（CEL＝5．0）＋生成項(xiàng)修正；

solver3：LU-SGS方法（CEL＝5．0）＋耗散項(xiàng)點(diǎn)隱式式計(jì)算＋生成項(xiàng)修正。

圖9 殘值隨迭代變化曲線（k-ωSST模型）Fig.9 Residual vs.number of iterations（k-ωSST）

圖10 翼型表面壓力系數(shù)分布曲線（k-ωSST模型）Fig.10 Distribution of pressure coefficient（k-ωSST）

從計(jì)算結(jié)果看出，只有solver3實(shí)現(xiàn)了較好的模擬，solver1和solver2都很快發(fā)散。這說明復(fù)雜的源項(xiàng)給k-ωSST模型控制方程組帶來了很大的剛性，必須結(jié)合采用多種手段處理才能得到理想的結(jié)果。另外，同樣是標(biāo)準(zhǔn)網(wǎng)格算例，k-ωSST模型收斂速度明顯較S-A模型慢，這也佐證了由于剛性較大而導(dǎo)致的局部時(shí)間步長進(jìn)一步減小的結(jié)論。

4 結(jié) 論

本文以S-A一方程模型和k-ωSST兩方程模型為例，研究了湍流模型方程源項(xiàng)對控制方程組產(chǎn)生的影響以及在實(shí)際應(yīng)用中的求解辦法。源項(xiàng)對原控制方程組雙曲性態(tài)的影響使得湍流變量在時(shí)間尺度和空間尺度的劇烈變化，帶來了嚴(yán)重的數(shù)值剛性問題。對此除了將源項(xiàng)和時(shí)間導(dǎo)數(shù)項(xiàng)耦合點(diǎn)隱式處理這一常規(guī)思路外，必須綜合考慮生成項(xiàng)和耗散項(xiàng)之間的平衡關(guān)系（特別是在邊界層中），防止迭代誤差導(dǎo)致過大或過小的源項(xiàng)值的擴(kuò)散——這需要從流動物理機(jī)制出發(fā)，對不同湍流模型源項(xiàng)做不同的分析處理。數(shù)值驗(yàn)證結(jié)果顯示無論是對源項(xiàng)剛度較小的S-A模型還是較大的k-ωSST模型，在有針對性的源項(xiàng)處理后都得到了良好的收斂效果。

［1］YAN C．Methods and applications of computational fluiddymanics［M］．Beijing：Beihang University Press，2006．（in Chinese）

閻超．計(jì)算流體力學(xué)方法及應(yīng)用［M］．北京：北京航空航天大學(xué)出版社，2006．

［2］LI X L，EU D X，MA Y W．Direct numerical simulation of hypersonic boundary layer transition over a blunt cone［J］．AIAA Journal，2008，46（11）：2899-2913．

［3］HUANG P G，COAKLEY T J．An implicit Navier-Stokes code for turbulent flow modeling［R］．AIAA 92-0547．

［4］BLAZEK J．Computational fluiddynamics：principles and applications［M］．London：ELSEVIER Science Ltd，2001．

［5］DU T，WU Z N．Mixed analytical／numerical method for low Reynolds-numberk-epsilon turbulence models［J］．AIAA Journal，2004，42：1140-1153．

［6］DU T，WU Z N．Mixed analytical／numerical method applied to the high Reynolds numberκ-εturbulence model［J］．Com puters &Eluids，2005，34（1）：97-119．

［7］DU T，WU Z N．Mixed analytical／numerical method for flow equations with a source term［J］．Computers&Eluids，2003，32（5）：659-690．

［8］DU T，YANG Y，WU Z N．Mixed analytical／numerical method and its application to problems of turbulent flow simulation［J］．Acta Aeronautica Et Astronautica Sinica，2007，27（2），198-203．（in Chinese）

杜濤，楊勇，吳子牛．混合解析／數(shù)值方法及其在湍流數(shù)值模擬上的應(yīng)用［J］．航空學(xué)報(bào)，2007，27（2）：198-203．

［9］SPALART P R，ALLMARAS S R．A one-equation turbulence model for aerodynamic flows［R］．AIAA 920439．

［10］COOK P H，MCDONALD M A，EIRMIN M C．Aerofoil RAE 2822 pressuredistributions，and boundary layer and wake measurements［R］．Experimental Data Base for Computer Program Assessment，AGARD AR 138，1979．

［11］MENTER E R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32（8）：1598-1605．

［12］MERCI B，STEELANT J，VIERENDEELS J．Computational treatment of source terms in two-equation turbulence models［J］．AIAA Journal，2000，38（11）：2085-2093．

Effects of numerical method of source terms in turbulence models

WANG Xiangyu，LI Dong
（National Key Laboratory of Science and Tecnology on Aerod ynamic Design and Research，Northwestern Polytecnical University，Xi’an 710072，China）

Eddy viscosity models have been widely used in many engineering applications of turbulence．Usually there are additional source terms on the right side of the turbulent governing equations，which may play an important role in the computation of Navier-Stokes equations．Stiffness that results from source terms is analyzed indetails and the significant balance between production term anddestruction term isdiscussed，together with point implicit approach，an universal way to solve equations with source terms．Afterwards both S-A turbulence model with weak stiffness andk-ωSST turbulence model with strong one are taken as examples to show the specific treatments concerning about thedifferent source terms between these two．The paper ends with numerical tests of RAE 2822 airfoil transonic flow，whichdemonstrates that the methods above can overcome the stiffness from the source terms and improve the stability of the computation effectively．

Navier Stokes equations；source term；stiffness；S-A turbulence model；k-ωSST turbulence model；computer simulation；airfoils

V211.3

Adoi：10.7638／kqdlxxb-2012.0224

0258-1825（2014）06-0868-06

2013-01-05；

2013-05-13

國家自然科學(xué)基金（11072200）

王翔宇（1989-），男，河北保定人，博士研究生，研究方向：計(jì)算流體力學(xué)．E-mail:xiangyu.wang＠m(xù)ail.nwpu.edu.cn

王翔宇，李棟．湍流模型源項(xiàng)的數(shù)值處理方法對計(jì)算結(jié)果的影響研究［J］．空氣動力學(xué)學(xué)報(bào)，2014，32（6）：868-873．

10．7638／kqdlxxb-2012．0224 WANG X Y，LI D．Effects of numerical method of source terms in turbulence models［J］．ACTA Aerodynamica Sinica，2014，32（6）：868-873．