網(wǎng)格和物理模型在計算流體動力學(xué)中的應(yīng)用研究

張亞勃，張東輝

（1.中國中原對外工程有限公司，北京100191；2.中國原子能科學(xué)研究院，北京102413）

1 CFD數(shù)值分析的基礎(chǔ)：網(wǎng)格生成技術(shù)

作為現(xiàn)代研究手段的數(shù)值分析方法，無論研究領(lǐng)域如何，研究區(qū)域的離散，即網(wǎng)格的生成是必不可少的一部分，根據(jù)研究問題的復(fù)雜程度的不同，這部分的工作量幾乎與進行求解計算的工作量相當(dāng)。因此合理有效的網(wǎng)格生成技術(shù)非常重要。由于工作量很大，在網(wǎng)格生成方法上的選擇主要以能生成合理的網(wǎng)格為原則，兼顧方法的簡便性，同一控制體中也可采用不同區(qū)域靈活選取不同網(wǎng)格生成方法。網(wǎng)格生成方法大體可分為結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。結(jié)構(gòu)化網(wǎng)格中，任意一個節(jié)點的位置可以通過一定的法則予以命名，當(dāng)研究區(qū)域復(fù)雜時，不同區(qū)網(wǎng)格按照一定規(guī)律聯(lián)接，形成塊結(jié)構(gòu)化網(wǎng)格，將控制體分解成不同的子域，每一子域中結(jié)構(gòu)化形式可以不同，子域之間可以互不重疊也可以有部分重疊。非結(jié)構(gòu)化網(wǎng)格中，節(jié)點位置無法用一個固定的法則予以有序的命名，非結(jié)構(gòu)化網(wǎng)格由于對不規(guī)則區(qū)域的特別適應(yīng)性而自20世紀80年代以來得到迅速發(fā)展，這種網(wǎng)格除了每一個單元及節(jié)點的幾何信息必須儲存外，與相鄰單元的編號等信息也必須作為聯(lián)系關(guān)系信息存儲起來，這使得非結(jié)構(gòu)化網(wǎng)格的信息量比較大。目前主流的CFD分析軟件都兼容兩種網(wǎng)格形式，并提供了各種各樣的網(wǎng)格生產(chǎn)方法，在同一問題中還可同時應(yīng)用這兩種不同形式的網(wǎng)格混合使用，達到特殊形狀模擬和提高計算效率的目的。

反應(yīng)堆熱工水力CFD模擬與傳統(tǒng)的主要分析復(fù)雜外流場流動的航空空氣動力學(xué)、內(nèi)燃機械等不同的是：反應(yīng)堆熱工水力CFD分析的主要對象是被反應(yīng)堆冷卻劑所包圍的反應(yīng)堆內(nèi)部構(gòu)件，如堆芯組件、控制棒組件、主泵等，這些特殊結(jié)構(gòu)的特點是設(shè)備全部或部分浸泡在冷卻劑中，而同時設(shè)備內(nèi)部也存在流體的流動，即在分析過程中既要考慮外流場的復(fù)雜流動現(xiàn)象，又要考慮內(nèi)部流場的影響；同時在計算中普遍要耦合熱效應(yīng)。這給反應(yīng)堆設(shè)備CFD網(wǎng)格模型的劃分帶來新的問題和挑戰(zhàn)。近年來，基于各種不同類型的結(jié)構(gòu)／結(jié)構(gòu)雜交網(wǎng)格技術(shù)的開發(fā)，隨著對復(fù)雜形狀問題適應(yīng)能力越來越迫切，一種新型的網(wǎng)格技術(shù)——多面體網(wǎng)格技術(shù)逐漸成為成熟的網(wǎng)格技術(shù)。在復(fù)雜網(wǎng)格生成與計算效率上，多面體網(wǎng)格提供了很好的平衡，相對于結(jié)構(gòu)網(wǎng)格，多面體網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格一樣簡單、快捷的構(gòu)建，而計算效率明顯高于非結(jié)構(gòu)化網(wǎng)格。對于應(yīng)用于內(nèi)流場形狀十分復(fù)雜的結(jié)構(gòu)，顯示出了較大的適用性。目前在反應(yīng)堆熱工流體設(shè)計中，對堆芯組件內(nèi)通道的模擬已經(jīng)廣泛地被應(yīng)用（如圖1所示）。

圖1 某快堆燃料組件多面體網(wǎng)格模型Fig.1 Polyhedral meshed model for fast reactor assembly

2 數(shù)值分析的離散格式：有限容積法

對流-擴散問題數(shù)值計算的精度與效率，不僅取決于所生成網(wǎng)格的質(zhì)量，同時也取決于對流項的離散格式（或差分格式）。在CFD研究領(lǐng)域里，對流項離散格式的構(gòu)造一直是具有挑戰(zhàn)性的研究熱點之一，長期以來認為難以構(gòu)造一種很好的模擬所有對流-擴散問題的對流項離散格式，因為在構(gòu)造對流項離散格式時不僅要考慮對流項遷移特性，還要保證較高的計算精度，這本身是互相對立的。低階格式如一階迎風(fēng)（FUD）格式是絕對穩(wěn)定的，能很好地體現(xiàn)對流項的遷移特性，但有非常嚴重的假擴散。中心差分（CD）格式能很好地消除假擴散，但它是條件穩(wěn)定的，當(dāng)網(wǎng)格貝克萊數(shù)或網(wǎng)格雷諾數(shù)超過臨界值時容易導(dǎo)致迭代發(fā)散。考慮到反應(yīng)堆組件CFD分析具有特殊形狀和邊界條件，且需要耦合熱效應(yīng)分析，因此有必要選用一種高分辨率的格式。

就常用的數(shù)值計算方法而言，大致可分為三類：有限差分法、有限容積法和有限元法。在CFD發(fā)展初期，由于所研究問題的幾何外形簡單，容易生成規(guī)則的結(jié)構(gòu)網(wǎng)格，可直接利用差分來近似微分方程的偏導(dǎo)數(shù)，因此有限差分方法得到很廣泛的應(yīng)用，隨著CFD面對計算的外形越來越復(fù)雜，傳統(tǒng)網(wǎng)格已經(jīng)不能適應(yīng)計算的要求，于是CFD工作者提出了分塊結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)化網(wǎng)格以及上文提到的多面體網(wǎng)格等，在這些網(wǎng)格上，傳統(tǒng)的有限差分法已經(jīng)不能再適用，為此相繼發(fā)展了有限容積法和有限元法。

有限元法是基于變分原理提出的計算方法，在計算固體力學(xué)中得到廣泛應(yīng)用，CFD方面起步較晚，其主要原因是流體力學(xué)的控制方程較固體力學(xué)更為復(fù)雜，且方程性質(zhì)隨流場特性（如雷諾數(shù)、馬赫數(shù)等）的不同呈現(xiàn)不同的性狀，變分原理不好應(yīng)用。有限元法求解難度和工作量較大，因此，目前來說有限元法的結(jié)算效率要低于有限容積法。

有限容積法通過對原方程（Euler和NS方程）在控制體單元內(nèi)直接積分，并通過體積分到面積分的轉(zhuǎn)化，將原方程的求解轉(zhuǎn)化為表明通量積分；黏性項和對流項也可轉(zhuǎn)化為通量面積分，從而求得單元內(nèi)物理量平均值。由于對控制體性狀沒有特殊的要求，此方法特別適合于復(fù)雜形狀流場的計算。因此，是進行反應(yīng)堆內(nèi)部流場計算的最佳選擇。

以直角坐標系中的對流項為例，它的有限體積離散形式為：

可由一階迎風(fēng)格式計算，對于對流占優(yōu)的流動，迎風(fēng)格式要優(yōu)于中心差分格式，此時：

3 湍流模型：特殊形狀控制體內(nèi)流場分析處理方法

工程角度研究湍流的基本目的是精確預(yù)測與控制湍流。隨著計算機技術(shù)的迅速發(fā)展，湍流的數(shù)值模擬精度幾乎可以與理論研究和實驗研究相比較，成為工程技術(shù)領(lǐng)域主要的分析手段。從求解方程上，湍流模型可以分為兩類：一種是以NS方程為求解對象，一種則以分子動力學(xué)理論研究湍流流動。前者主要分為直接模擬（direct numerical simulation，DNS）、大渦模擬（large eddy simulation，LES）、雷諾平均法（Reynolds averaged Navier-Stockes，RANS）和基于概率方法的PDF（probability density function）。DNS是對非穩(wěn)態(tài)NS方程直接求解，可以分辨所有尺度的渦，RANS求解雷諾平均方程，雷諾平均方程是不封閉的，必須引入雷諾應(yīng)力的封閉模型才能進行求解。LES是介于DNS和RANS之間的數(shù)值模擬方法，對大尺度渦進行直接求解，小渦對大渦的影響采用經(jīng)驗?zāi)Ｐ徒铺幚怼DF則應(yīng)用于湍流化學(xué)模擬反應(yīng)方面較多。目前工程應(yīng)用為目的的CFD分析大部分仍以雷諾平均法為主，雖然雷諾平均法需要引入較多經(jīng)驗參數(shù)修正，且分析精度相對低，但其優(yōu)點是經(jīng)過數(shù)十年的研究和應(yīng)用，雷諾平均法已經(jīng)發(fā)展得相當(dāng)成熟，經(jīng)過修整可以精確的模擬各種特性的流場，且計算效率高，能適應(yīng)復(fù)雜流場，巨大網(wǎng)格數(shù)量的問題分析。因此，組件及冷卻劑熱工流體動力學(xué)設(shè)計和分析還是以雷諾平均法為主，而對于形體復(fù)雜的反應(yīng)堆系統(tǒng)來說，進行雷諾平均法模擬是比較現(xiàn)實的。下面介紹幾種實踐中被認為較為適合的特殊流道湍流模型，對于更多的湍流模型可見圖2，詳細介紹可參考文獻［2］。

3.1 矩形截面通道內(nèi)湍流的大渦模擬

大渦模擬介于雷諾平均法模擬和DNS直接模擬之間，它既能反映控制體內(nèi)大部分渦的運動特性，又簡化了數(shù)量眾多且影響較小的小尺度渦的直接模擬，提高計算效率，同時其計算精度要較雷諾平均法模擬高。近年來的發(fā)展，該模型能夠?qū)Σ灰?guī)則形體進行較高效率的計算。在不規(guī)則管道或流場中，由于湍流的各向異性引起垂直于主流方向的二次流動，使得湍流和換熱更加復(fù)雜。對于矩形管道內(nèi)湍流的模擬，各向同性的渦黏模型不能預(yù)測湍流引起的二次流動，代數(shù)應(yīng)力模型、非線性k-ε模型和各向異性低雷諾數(shù)模型經(jīng)常被使用，但對二次流動的表達取決于經(jīng)驗參數(shù)，因此模擬出的湍流量與試驗值不能完全一致。為更好地理解和揭示湍流機理和二次流動產(chǎn)生原因，直接模擬是一個重要的工具。許多研究者采用大渦模擬技術(shù)研究矩形通道內(nèi)的流動和換熱情況，這些研究成功地預(yù)測了截面內(nèi)二次流動以及對平均流動和湍流統(tǒng)計的影響，且與試驗是一致的［3］。

例如某矩形管道高為H，長為L，x為主流方向，y和z分別為展向，u、v、w分別表示各方向的速度分量，考慮重力影響，其大渦模擬物理模型描述如下：

式中字母上方的“—”表示過濾量，無量綱溫度的定義為：

根據(jù)理論抽樣，本文選擇在參與社區(qū)治理之初規(guī)制、規(guī)范和認知合法性均缺乏的LL老年社會工作服務(wù)中心(以下簡稱LL)作為案例。作為純粹的“草根”(自下而上發(fā)起成立的)社會組織，LL成立于2006年，是一家致力于推動中國社區(qū)居家養(yǎng)老服務(wù)的民間機構(gòu)。從志愿團隊身份到企業(yè)法人身份，再到合法社會組織身份，并成為行業(yè)標桿組織，其合法化路徑與行動策略顯示出社會組織自下而上獲得合法性的完整過程，較好地排除了組織資源、政治聯(lián)系等其他因素對社會組織合法化進程的影響。

雷諾數(shù)和格拉曉夫數(shù)為：

以上式中τij和qj分別表示亞格子湍流應(yīng)力和熱流密度，大渦模擬中，亞格子模型對大渦模擬具有重要影響，若采用動態(tài)渦黏模型表達則亞格子應(yīng)力表示為：

以下為某對流實驗回路的模擬，裝置的主要結(jié)構(gòu)由右側(cè)高溫段、左側(cè)低溫段、下端加熱段和上端冷卻段組成。回路圓管直徑60.3mm，壁厚5mm，其結(jié)構(gòu)參數(shù)如圖3所示。

圖3 高溫實驗回路及其結(jié)構(gòu)Fig.3 High temperature test loop and its diagrammatic sketch

管內(nèi)流體為液態(tài)鋰鉛，其物性參數(shù)如下：

密度／（kg／m3）：

比熱容／（J／kg·K）：

動力黏度／（Pa·s）：

利用CFD程序?qū)Ω邷責(zé)釋α鲗嶒灮芈愤M行計算，計算結(jié)果見表1。回路上升段溫度沿高度方向上的分布如圖4所示。

表1 計算與實驗數(shù)據(jù)的對比Table 1 Comparison of calculated and experimental data

圖4 回路上升段溫度沿高度方向上的分布Fig.4 Temperature distribution along height direction

3.2 截面通道內(nèi)存在孤島湍流的雷諾平均模擬

雷諾平均法模型是在各向同性的前提下，將速度脈動的二階關(guān)聯(lián)量表示成平均速度梯度與湍流黏性系數(shù)的乘積，即：

推廣到三維問題，若用笛卡爾張量表示則有：

湍流模型或湍流封閉的任務(wù)就是給出計算湍流黏性系數(shù)的表達式或其輸運方程的方法。根據(jù)建立微分方程的數(shù)目，可以分為零方程模型（代數(shù)方程模型）、一方程模型和雙方程模型。由于該模型是基于各向同性的前提下，對于特殊形狀通道的二次流計算需要對模型本身進行修正，以達到模擬二次流對主流影響的效果，如修正的k-ε模型，低雷諾數(shù)k-ε模型等。雷諾平均法模型計算效率非常高，適合計算復(fù)雜形狀、網(wǎng)格數(shù)量大的問題。

在快堆液體懸浮式非能動停堆裝置的設(shè)計中，懸浮式控制棒懸浮于液鈉中間，屬于流固耦合計算的孤島問題。所謂孤島的處理方法，是在采用整場離散、整場求解時，對于與計算邊界相鄰接的固體區(qū)域可以采用把固體區(qū)作為黏性無限大的流體區(qū)的方法來保證固體區(qū)中速度為零。當(dāng)固體區(qū)位于流場中間而不與計算邊界相連時采用的處理方法。當(dāng)然該處理手段并不只限于處理孤立物體，而是一種通用性很強的處理流固耦合問題的方法。該方法具體實施步驟如下：①在每一層次的迭代計算前令障礙物區(qū)的速度u和為零，這樣可保證對障礙物附近的節(jié)點速度起到滯止的作用；②在求解速度的方程時令障礙物區(qū)的離散方程的主對角線上的系數(shù)為一個極大值，如1 024，這樣算出來的速度估計值u和v等于零；③再令障礙物區(qū)的速度修正公式里的系數(shù)d，d取一個極小值，如10～24，這樣速度修正值也為零。關(guān)于障礙物絕熱的邊界條件，采用的是使障礙物和其鄰近的流體溫差等于零的方法。即在每一層次迭代計算前，令障礙物邊界的節(jié)點上的溫度值等于其鄰近流體的溫度。至于上下底邊的絕熱邊界條件，按附加源項法處理時，其附加的常數(shù)源項為零。

同時，對于形狀特殊，流道狹窄的堆芯控制棒等裝置，在湍流模型的選取上應(yīng)選用二次流修正較好，附加壁面影響的模型。徐建軍等人采用染色法和數(shù)值計算相結(jié)合的方法研究冷態(tài)下矩形窄縫流道角部附近流場的分布，研究表明，相比較其他湍流模型，雷諾應(yīng)力模型中的SSG模型能夠模擬出矩形通道中湍流的各向異性，與可視化實驗結(jié)果定性符合較好。標準k-ε模型在低雷諾數(shù)時符合較好，但隨著雷諾數(shù)的增加，預(yù)測偏差增大。在較寬的雷諾數(shù)范圍內(nèi)，SSG模型與試驗數(shù)據(jù)擬合的公式符合較好，尤其是隨著Re的增加，偏差減小，在工程應(yīng)用中，控制棒組件中的Re比較高，選用SSG模型是合適的，見圖5。

圖5 試驗數(shù)據(jù)與CFD計算結(jié)果比較Fig.5 Comparison of calculated and experimental data

4 發(fā)展趨勢：兩相流、多尺度耦合和并行計算的應(yīng)用

兩相流動和沸騰傳熱的計算分析是一項復(fù)雜的工作，因為通過氣液界面的流動特性和傳熱特性變化較大，特別是密度的巨大變化和熱導(dǎo)率的改變：與運用一維集總參數(shù)法進行兩相流模擬的系統(tǒng)程序相比，CFD方法是基于求解空間兩相流場內(nèi)各相的質(zhì)量守恒、動量守恒和能量守恒方程。兩相流的模擬還無法達到單相模擬所能達到的精度。總之，模擬兩相流動比單相流動CPD模擬更富有挑戰(zhàn)性，但其潛力是巨大的，目前兩相流動CPD模擬方法正在緊張而有序的研究發(fā)展中。

早在20世紀80年代Hirt等人提出了所謂流體域VOF（Volume Of Fluid）方法直接進行兩組分、兩相界向邊界的數(shù)值模擬，但沒有考慮相變熱。VOF方法已經(jīng)成為幾乎所有模擬兩相流動商用CFD程序的基礎(chǔ)。近年來，隨著多相流物理模型特別是兩流體模型的發(fā)展，多相流動的數(shù)值模擬方法得到一些應(yīng)用，開發(fā)了許多可用于模擬多相流動的CFD軟件：CFX，lquent.FLOW3D，STAR—CD，Phoenix等。多相流的描述主要考慮總體的性能：壓降、體積份額和傳熱率等。目前，根據(jù)主要模擬方法的原理，多相流動模擬方法主要包括直接數(shù)值模擬方法DNS（direct numerical simulation）、均相模型、歐拉—拉格朗日方法和兩流體模型等。

盡管在工程和研究方面取得了一些進展，但對以兩流體模型為基礎(chǔ)的多維兩相流的數(shù)值模擬仍處于發(fā)展階段，并且流入流出邊界很可能引起數(shù)值上困難和收斂問題。工程應(yīng)用主要需考慮魯棒性、CPU時間效率和準確性。并且多相流動的模擬不像單相流那樣，僅提高網(wǎng)格的質(zhì)量就能改善計算結(jié)果的準確性，物理模型的影響更為重要。

經(jīng)過近半個世紀的研究，特別是近十年來的發(fā)展，CFD方法得到大力的發(fā)展和廣泛的應(yīng)用而且針對具體問題，如核工業(yè)中的問題必將推動CFD方法的進一步發(fā)展和完善，CFD方法使計算分析由一維向三維發(fā)展，并且三維CFD計算往往與一維系統(tǒng)程序計算進行耦合，即所謂多尺度耦合。現(xiàn)有許多項目就致力于類似的工作：例如，愛達荷州工程和環(huán)境國家實驗室已與FLUENT軟件公司聯(lián)手將RELAP5—3D和FLUENT進行多尺度耦合。歐洲的ASTAR項目的目的就是在于將先進的CFD工具與系統(tǒng)程序（CATHARE或ATHLET）進行耦合，德國也正將CFX模擬和系統(tǒng)程序COCOSYS的計算進行耦合。近來就有學(xué)者使用多尺度耦合方法對反應(yīng)堆堆芯進行研究。D.P.Weber等在2003年對微型堆芯（mini core）進行了耦合計算，使用了5 000 000CFD網(wǎng)格單元和150 000中子物理學(xué)區(qū)域。初步計算表明，如果采用萬億次級別的計算機，原型PWR堆芯的穩(wěn)態(tài)計算能在一天內(nèi)完成。但系統(tǒng)程序和CFD的耦合，僅限于穩(wěn)態(tài)工況，系統(tǒng)程序為詳細的CFD分析提供邊界條件。高速和大存儲量的計算機的發(fā)展使運用CFD方法可得到許多流動問題的數(shù)值解。

在核工業(yè)系統(tǒng)中，許多問題較為復(fù)雜，需要分析的對象幾何結(jié)構(gòu)尺寸較大，采用CFD方法進行研究時往往受到計算機計算能力的限制，因此必須采用并行計算方法。目前并行計算面臨兩個主要問題，一是如何將龐大的數(shù)據(jù)分配到不同的中央處理器上進行計算；二是并行計算方法的魯棒性較差，Weber等人1999年用STARCD對單相、穩(wěn)態(tài)情況下1／8個PWR堆芯內(nèi)的溫度分布進行了計算，以驗證汽泡核心可能出現(xiàn)的位置，整個區(qū)域劃分了2.5億個網(wǎng)格單元，采用計算機群進行并行計算，如果要對該系統(tǒng)進行瞬態(tài)計算，則只允許劃分數(shù)百萬的網(wǎng)格。

5 不足與挑戰(zhàn)

CFD方法應(yīng)用于反應(yīng)堆工程設(shè)計的過程中取得了很大的成就，但它畢竟還是一門年輕的學(xué)科，存在著很多困惑。

首先是可信度問題。對于絕大多數(shù)流動問題，還沒有任何一個CFD方法可以給出定性的正確結(jié)論，或確定誤差范圍定量結(jié)果。一般采用同實驗結(jié)果相比較的方法來驗證CFD的可信度問題，但一方面，試驗結(jié)果同樣存在可信度問題，另一方面，實驗結(jié)果往往很難給出精細的流動結(jié)構(gòu)，而對細節(jié)的深入分析是CFD結(jié)果的重要信息。

其次，CFD理論的數(shù)學(xué)基礎(chǔ)尚不健全。目前非線性偏微分方程組的數(shù)值解法尚未成熟，對三維流體力學(xué)建立的控制方程在進行計算處理時必須是經(jīng)過簡化的，Euler和N-S方程本身都是非線性方程組，針對標量方程構(gòu)造的CFD計算格式通過特征分裂等方法擴展到這些非線性方程組后，其原來的一些特性就消失了。CFD的其他一些物理模型如湍流模型、化學(xué)反應(yīng)模型、多相流模型等都遠未發(fā)展完善，其計算精度和可靠性都存在目前無法解決的缺陷，解決這些問題可能需要很長時間的探索。

6 結(jié)語

越來越多的研究表明，今后CFD方法將由反應(yīng)堆設(shè)計安全審查部門要求應(yīng)用于安全評審。目前在開發(fā)新一代的反應(yīng)堆時已經(jīng)使用該研究方法，在第四代壓水堆概念設(shè)計時，將應(yīng)用CFD程序評估完整的壓力容器特性，綜合分析壓力容器的熱量和氣體傳輸以及其對系統(tǒng)壓力變化的反應(yīng)性。隨著現(xiàn)代計算機技術(shù)特別是超級計算機和并行計算集群等技術(shù)的發(fā)展，CFD方法逐漸開始擺脫計算機硬件的限制，其應(yīng)用范圍將大大擴展和豐富，成為反應(yīng)堆工程設(shè)計與熱工水力安全分析的主要手段。

［1］ 王福軍.計算流體動力學(xué)分析［M］.北京：清華大學(xué)出版社，2004.

［2］ 林誠格，等.非能動安全先進核電廠AP1000［M］.北京：原子能出版社，2008.

［3］ Patankar，S.V.（1975），Numerical Prediction of Three-Dimensional Flows，Studies in Convection Theory，Measurement and Applications，Vol.1Academic，New York.

［4］ Fluent Inc.FLUENT User's Guide.Fluent Inc.2003.