網格和物理模型在計算流體動力學中的應用研究

張亞勃，張東輝

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

1 CFD數值分析的基礎：網格生成技術

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

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

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

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

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

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

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

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

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

可由一階迎風格式計算，對于對流占優的流動，迎風格式要優于中心差分格式，此時：

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

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

3.1 矩形截面通道內湍流的大渦模擬

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

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

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

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

雷諾數和格拉曉夫數為：

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

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

圖3 高溫實驗回路及其結構Fig.3 High temperature test loop and its diagrammatic sketch

管內流體為液態鋰鉛，其物性參數如下：

密度／（kg／m3）：

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

動力黏度／（Pa·s）：

利用CFD程序對高溫熱對流實驗回路進行計算，計算結果見表1。回路上升段溫度沿高度方向上的分布如圖4所示。

表1 計算與實驗數據的對比Table 1 Comparison of calculated and experimental data

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

3.2 截面通道內存在孤島湍流的雷諾平均模擬

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

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

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

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

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

圖5 試驗數據與CFD計算結果比較Fig.5 Comparison of calculated and experimental data

4 發展趨勢：兩相流、多尺度耦合和并行計算的應用

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

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

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

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

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

5 不足與挑戰

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

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

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

6 結語

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

［1］ 王福軍.計算流體動力學分析［M］.北京：清華大學出版社，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.