符棟良

（上海市特種設(shè)備監(jiān)督檢驗(yàn)技術(shù)研究院）

工程上有許多裝置都很龐大，有的還處于高溫、高壓甚至有害的環(huán)境下工作，因此要進(jìn)行有關(guān)的實(shí)驗(yàn)往往很困難，甚至是不太現(xiàn)實(shí)的。于是我們就用數(shù)值模擬，即用數(shù)學(xué)的方法模擬實(shí)際問題。隨著計(jì)算機(jī)存儲(chǔ)容量的不斷增加，計(jì)算速度的不斷提高，計(jì)算方法的不斷完善，數(shù)值模擬得到了迅速的發(fā)展。然而對(duì)大規(guī)模問題進(jìn)行數(shù)值模擬往往還需要消耗大量的計(jì)算機(jī)資源，所以有必要用多機(jī)來進(jìn)行并行求解以提高計(jì)算效率。

氣固兩相流動(dòng)廣泛存在于各種分離和燃燒裝置中，正確預(yù)報(bào)其中兩相速度、顆粒濃度及兩相湍動(dòng)能分布，對(duì)了解流動(dòng)、混合及反應(yīng)機(jī)理，進(jìn)行優(yōu)化設(shè)計(jì)以提高性能，都很重要?；鹆Πl(fā)電廠燃燒器中空氣和煤粉混合物的流動(dòng)就是典型的多相流動(dòng)工況[1]。其計(jì)算規(guī)模很大，加之其計(jì)算區(qū)域比較復(fù)雜，所以對(duì)并行計(jì)算更是有迫切的需求。

1 區(qū)域分解法

區(qū)域分解法[2]是把計(jì)算區(qū)域分解為若干個(gè)子區(qū)域，子區(qū)域應(yīng)盡可能規(guī)則，從而將原問題的求解轉(zhuǎn)化為在子域上的求解。它把大問題化為若干個(gè)小問題，從而縮小了計(jì)算規(guī)模。這對(duì)于克服計(jì)算機(jī)容量和速度的限制具有極大的意義。本論文正是基于這一想法，決定采用區(qū)域分解算法，實(shí)現(xiàn)對(duì)燃燒器一次風(fēng)噴咀和彎頭組件內(nèi)氣固兩相流動(dòng)問題的分區(qū)并行求解。根據(jù)區(qū)域分解理論，各子區(qū)域劃分好以后，關(guān)鍵的問題是子區(qū)域之間交界面上的信息如何進(jìn)行快速有效的傳遞，保證計(jì)算的收斂性。這也就是所謂的交界面技術(shù) （interface method）。

交界面技術(shù)包括兩大類，即非重疊區(qū)域耦合和重疊區(qū)域耦合。前者具有簡(jiǎn)單方便的優(yōu)點(diǎn)，但是應(yīng)用范圍不夠廣。后者應(yīng)用范圍廣，但是一般要求耦合區(qū)域具有相同的網(wǎng)格系統(tǒng)且不得劃分更多的子區(qū)域，這對(duì)于復(fù)雜區(qū)域來說顯得不夠靈活。本文所采用的是非重疊區(qū)域分解算法。根據(jù)區(qū)域分解的思想和算法，將計(jì)算區(qū)域進(jìn)行合理的劃分，可以實(shí)現(xiàn)分區(qū)計(jì)算。若進(jìn)一步結(jié)合并行技術(shù)，則可利用現(xiàn)有的微機(jī)資源，在微機(jī)網(wǎng)絡(luò)上實(shí)現(xiàn)較大規(guī)模問題的分區(qū)并行計(jì)算

2 并行計(jì)算環(huán)境

2.1 MPI介紹

本文采用的網(wǎng)絡(luò)并行環(huán)境是MPI(message passing interface)，MPI是目前最重要的實(shí)現(xiàn)消息傳遞的并行編程環(huán)境之一，它具有移植性好、功能強(qiáng)大、效率高等多種優(yōu)點(diǎn)[3]。具體地講，MPI就是提供了并行庫。這樣，原來的串行編譯器能夠繼續(xù)使用，編程者只需要在原來的串行程序中加入對(duì)并行庫的調(diào)用，就可以實(shí)現(xiàn)并行程序設(shè)計(jì)。

一個(gè)典型的基于MPI并行環(huán)境的程序框架結(jié)構(gòu)可以用圖1來表示。

圖1 MPI程序的框架結(jié)構(gòu)

2.2 MPI并行程序的主從模式

MPI的兩種最基本的并行程序設(shè)計(jì)模式是主從模式和對(duì)等模式。本文所選用的是主從模式，這種計(jì)算模式有一個(gè)控制程序稱為主進(jìn)程，負(fù)責(zé)進(jìn)程的生成、發(fā)送、收集數(shù)據(jù)并顯示計(jì)算結(jié)果，其余的從進(jìn)程執(zhí)行實(shí)際計(jì)算，如圖2所示。

圖2 并行計(jì)算主從模式的流程

可以說，服務(wù)器提供服務(wù)的軟件或進(jìn)程，而客戶機(jī)則接受服務(wù)的軟件或進(jìn)程。服務(wù)器和客戶機(jī)是針對(duì)網(wǎng)絡(luò)的通信進(jìn)程而言的，不是針對(duì)網(wǎng)絡(luò)客戶或個(gè)人計(jì)算機(jī)的。區(qū)分服務(wù)器和客戶機(jī)主要取決于通信進(jìn)程雙方的地位。同一臺(tái)個(gè)人計(jì)算機(jī)，用戶設(shè)置不同的進(jìn)程，那么用戶可能既是某個(gè)進(jìn)程的客戶機(jī)，又是另一個(gè)進(jìn)程的服務(wù)器。

主從模式具有以下兩個(gè)優(yōu)點(diǎn)。

（1）主從模型可以應(yīng)付網(wǎng)絡(luò)中資源、運(yùn)算能力和信息的不均等現(xiàn)象。

主從模式的最重要的特點(diǎn)是非對(duì)等的相互作用，即主進(jìn)程和從進(jìn)程處于不平等的地位。服務(wù)器擁有客戶機(jī)所不具備的硬件、軟件資源和運(yùn)算能力，服務(wù)器提供服務(wù)，客戶機(jī)請(qǐng)求服務(wù)。同時(shí)，網(wǎng)絡(luò)中還存在人為的不均等現(xiàn)象。這種不均等主要是指信息而言，例如為節(jié)約網(wǎng)絡(luò)資源沒有必要在每臺(tái)機(jī)上都保存一個(gè)拷貝，特別是在科學(xué)計(jì)算中存在的動(dòng)態(tài)變化信息，如果每臺(tái)主機(jī)都保存一份拷貝，那么為保持所有拷貝的一致性，要進(jìn)行的動(dòng)態(tài)修改和刷新工作，其開銷是非常巨大的。

（2）主從模式可以使數(shù)據(jù)交換保持同步。

網(wǎng)絡(luò)通信不同于單個(gè)計(jì)算機(jī)之間的通信的一個(gè)特點(diǎn)是，網(wǎng)絡(luò)通信是完全異步的，相互通信之間既不是父子關(guān)系，也不共享內(nèi)存緩沖區(qū)。因此需要一種機(jī)制，為準(zhǔn)備通信的進(jìn)程之間建立聯(lián)系，使二者的數(shù)據(jù)交換保持同步。而主從模型解決了上述問題。在主從模式中，每次通信均由隨機(jī)啟動(dòng)的從進(jìn)程發(fā)出，主進(jìn)程一直處于等待狀態(tài)。這樣可以保證主進(jìn)程隨時(shí)對(duì)從進(jìn)程請(qǐng)求做出響應(yīng)。

3 燃燒器的物理模型

本文的研究對(duì)象是貴州某電廠300 MW機(jī)組中的鍋爐用燃燒器，即1025 t/h亞臨界壓力無煙煤鍋爐中使用的一種新型燃燒器。

這是一個(gè)四角切圓直流燃燒器，和一般燃燒器相比，新型燃燒器在燃燒器一次風(fēng)噴咀靠近出口處增加了一個(gè)菱形的障礙塊，其目的是使燃燒器出口煤粉較均勻。攜帶著煤粉的一次風(fēng)從垂直圓管的下部進(jìn)入燃燒器，從水平方管的右側(cè)出口噴入爐膛。管道除了垂直圓管和水平方管外，還有兩個(gè)彎曲段和一個(gè)圓管變方管過渡段。第一個(gè)彎曲段為45°，位于Y-Z平面上；第二個(gè)彎曲段為90°，是在與垂直平面和水平面都成45°的平面上。圓管直徑為426 mm，方管尺寸為寬580 mm、高406 mm。具體結(jié)構(gòu)如圖3所示。

圖3 帶有菱形體的一次風(fēng)噴嘴及彎頭組件

4 計(jì)算網(wǎng)格生成

燃燒器組件網(wǎng)格生成的整體思路是采用多塊技術(shù)，將物理域劃分成多個(gè)子域，在每個(gè)子域中利用貼體網(wǎng)格技術(shù)中的解橢圓形方程的方法生成結(jié)構(gòu)化網(wǎng)格。為了保證網(wǎng)格質(zhì)量，塊與塊交界面上的網(wǎng)格保持一致[4]。

在區(qū)域分解中，將整個(gè)計(jì)算區(qū)域分裂為六個(gè)子域Ωt，從入口到出口其名稱和網(wǎng)格規(guī)模依次為：第一部分 （包括垂直圓柱、90°彎管部分、圓管變方管的過渡段）65 mm×16 mm×16 mm、方管部分20 mm×16 mm×16 mm、菱形障礙塊的上部10 mm×16 mm×6 mm、菱形障礙塊的下部10 mm×16 mm×6 mm、出口部分的方管10 mm×16 mm×16 mm和出口后的大空間 100 mm×100 mm×100 mm。 ?Ωt(t=1，2，3，4，5，6)為其邊界，?Ω 為整個(gè)計(jì)算區(qū)域的邊界。圖4是燃燒器一次風(fēng)組件整體的網(wǎng)格分布圖。其中在數(shù)值模擬中用到的爐膛部分的網(wǎng)格圖未給出，僅給出燃燒器一次風(fēng)組件內(nèi)的網(wǎng)格圖。

圖4 燃燒器一次風(fēng)組件整體的網(wǎng)格分布圖

5 氣固兩相湍流數(shù)值解法

本文選擇單流體模型來模擬燃燒器一次風(fēng)噴咀及彎頭組件內(nèi)的氣固兩相流動(dòng)。

本文對(duì)氣固兩相流動(dòng)采用單向耦合[5-8]。顆粒相的求解是建立在氣相場(chǎng)得到粗收斂解的基礎(chǔ)上的。顆粒相的壓力近似取為氣相場(chǎng)的壓力；對(duì)于顆粒相的速度場(chǎng)和密度場(chǎng)，用冪函數(shù)差分法對(duì)固相動(dòng)量方程和連續(xù)性方程進(jìn)行離散，得到關(guān)于速度場(chǎng)和濃度場(chǎng)的代數(shù)方程組，然后耦合迭代求解，最終得到滿足顆粒相動(dòng)量方程和連續(xù)性方程的密度分布和速度分布。

6 并行方案

區(qū)域分解方法是一類適合于分區(qū)并行處理的偏微分方程數(shù)值求解的方法。它的基本思想就是，把物理問題的求解區(qū)域按照一定的原則 （如幾何形狀、物理特性等）分解成若干個(gè)相對(duì)簡(jiǎn)單的子區(qū)域，于是原問題的求解就轉(zhuǎn)化為在各子區(qū)域上的分別求解及各子區(qū)域之間的邊界耦合。

根據(jù)主從模式的特點(diǎn)，我們可以安排每個(gè)從進(jìn)程 （一臺(tái)客戶機(jī)運(yùn)行1至2個(gè)從進(jìn)程）負(fù)責(zé)一個(gè)子區(qū)域的計(jì)算。從進(jìn)程在迭代若干步后，首先要把一些相關(guān)的交界面信息傳給主進(jìn)程，然后主進(jìn)程會(huì)根據(jù)條件得到新一輪的交界面信息，最后再把這些新得到的信息傳給相關(guān)的從進(jìn)程。有關(guān)詳情可參照文獻(xiàn) [9]～文獻(xiàn) [15]。

本文對(duì)于燃燒器一次風(fēng)噴咀和彎頭組件內(nèi)的三維氣固兩相流動(dòng)問題使用5臺(tái)PC機(jī)并行求解，其中1臺(tái)PC機(jī)作為服務(wù)器運(yùn)行主進(jìn)程，另外4臺(tái)作為客戶機(jī)運(yùn)行從進(jìn)程。這5臺(tái)機(jī)器的軟硬件配置如下。

服務(wù)器： CPU，AMD 64×2 Dual； 內(nèi)存，1G；操作系統(tǒng)，Windows XP。

客戶機(jī) 1： CPU， AMD 64×2 Dual； 內(nèi)存， 1G；操作系統(tǒng)，Windows XP。

客戶機(jī) 2： CPU， AMD 64×2 Dual； 內(nèi)存， 1G；操作系統(tǒng)，Windows XP。

客戶機(jī) 3： CPU， AMD 64×2 Dual； 內(nèi)存， 1G；操作系統(tǒng)，Windows XP。

客戶機(jī) 4： CPU， AMD 64×2 Dual； 內(nèi)存， 1G；操作系統(tǒng)，Windows XP。

本文將整個(gè)燃燒器計(jì)算區(qū)域分成6個(gè)子塊，如圖4所示。每一塊即在計(jì)算過程中對(duì)應(yīng)于一個(gè)進(jìn)程。其中進(jìn)程1放在客戶機(jī)1上，進(jìn)程2、3放在客戶機(jī)2上，進(jìn)程4、5放在客戶機(jī)3上，進(jìn)程6放在客戶機(jī)4上。

這里的服務(wù)器主要負(fù)責(zé)交界面信息的耦合、數(shù)據(jù)轉(zhuǎn)發(fā)和判斷收斂等工作，而客戶機(jī)主要負(fù)責(zé)各計(jì)算區(qū)域邊界、入口和出口條件的初始化以及對(duì)一個(gè)（或多個(gè)）子區(qū)域的數(shù)值求解工作。圖5給出了服務(wù)器與客戶機(jī)的關(guān)系圖。其中S代表服務(wù)器，C1、C2、C3和C4代表4個(gè)客戶機(jī)，箭頭代表數(shù)據(jù)傳遞的方向。

圖5 服務(wù)器與客戶機(jī)關(guān)系

7 數(shù)值模擬結(jié)果及分析

以下對(duì)Re=105時(shí)燃燒器一次風(fēng)噴咀和彎頭組

件內(nèi)的氣固兩相流動(dòng)的速度、湍動(dòng)能、顆粒相濃度和并行效率等參數(shù)進(jìn)行分析。

7.1 流場(chǎng)分析

流體流經(jīng)菱形障礙塊后，同樣也會(huì)在障礙塊后方形成兩個(gè)尾渦區(qū)。需要指出，該尾渦區(qū)是有利于爐膛內(nèi)穩(wěn)定燃燒的。顆粒速度分布如圖6所示。

圖6 J平面的顆粒速度分布

7.2 整流現(xiàn)象

圖7為燃燒器一次風(fēng)噴咀水平管內(nèi)橫向截面的氣相速度矢量圖。圖7中y+為水平管軸線方向的無量綱化的刻度 （y+=y/din，y為水平管軸線方向的刻度，din為入口圓管的直徑）。

其具體代表的物理位置為：

y+=2.15—— 圓管第二次彎曲后的截面；

y+=3.00——圓管變方管過渡段的截面；

y+=4.65——水平管內(nèi)靠近菱形障礙塊前端的截面；

圖7 燃燒器內(nèi)橫向截面的氣體速度矢量圖

y+=6.15——水平管內(nèi)靠近菱形障礙塊后端的截面。

從圖7可以看到，當(dāng)流體從垂直管道經(jīng)過彎頭流入水平管道時(shí)，由于管道需要為流動(dòng)提供向心力而產(chǎn)生橫向壓差，從而導(dǎo)致流場(chǎng)在經(jīng)過彎頭后出現(xiàn)雙螺旋流現(xiàn)象，這就是所謂的 “彎頭效應(yīng)”。此外，類似于 “后向臺(tái)階效應(yīng)”，經(jīng)過彎頭的流動(dòng)會(huì)在內(nèi)側(cè)產(chǎn)生一定程度的分離并形成一個(gè)尺度較大的回流區(qū)，從而導(dǎo)致軸向速度沿橫向截面的分布很不對(duì)稱。然而，由于菱形障礙塊的存在，它有一定的整流作用，所以當(dāng)遇到障礙塊的時(shí)候，螺旋流現(xiàn)象得到消除，并且通過障礙塊后的流動(dòng)也變得上下比較對(duì)稱了。

7.3 湍動(dòng)能分析

圖8為燃燒器一次風(fēng)噴咀水平管內(nèi)障礙塊前后的Y-Z平面的湍動(dòng)能分布圖。

其中的湍動(dòng)能為無量綱湍動(dòng)能，其定義為k+=k/u2in，其中k代表實(shí)際的湍動(dòng)能，uin代表實(shí)際的氣體入口速度。x+為水平管徑向的無量綱化的刻度（x+=x/din，x為燃燒器水平管水平徑向的刻度，din為入口圓管的直徑）。

從圖8可以看出，由于障礙塊的出現(xiàn)，使得流場(chǎng)受到干擾，從而障礙塊前后的湍動(dòng)能有明顯的提高。

7.4 濃度場(chǎng)分析

圖8 燃燒器內(nèi)障礙塊前后Y-Z平面的湍動(dòng)能分布

圖9為燃燒器一次風(fēng)噴咀近出口處Y-Z平面內(nèi)的顆粒相濃度分布圖。圖9中的濃度值為無量綱值，其定義為nk+=nk/nkin，其中nk為實(shí)際顆粒相的濃度，nkin為實(shí)際顆粒相的入口濃度。x+為燃燒器一次風(fēng)噴咀水平管水平徑向的無量綱化的刻度 （x+=x/din，x為水平管水平徑向的刻度，din為入口圓管的直徑）。

從圖9我們明顯可以看出，由于受到障礙塊阻擋，煤粉在障礙塊前有堆積，所以濃度較高；由于該菱形障礙塊的布置位置離彎頭較遠(yuǎn)，所以不能起到濃淡分離的作用。

圖9 燃燒器一次風(fēng)噴咀近出口處Y-Z平面的顆粒相濃度分布

7.5 并行性能分析

下文從并行計(jì)算的并行執(zhí)行時(shí)間、加速比和并行效率上來進(jìn)行并行性能的分析。

其計(jì)算公式如下：

式中Sn——加速比，衡量并行機(jī)群有多快；

t1——順序執(zhí)行時(shí)間，即程序在單處理機(jī)上的運(yùn)行時(shí)間；

tn——并行執(zhí)行時(shí)間，即程序在并行機(jī)群上的運(yùn)行時(shí)間；

En——并行效率，衡量處理機(jī)的利用率；

n——處理器的數(shù)目。

表1列出了使用單機(jī)求解和多機(jī)并行求解該問題的并行性能參數(shù)。其中多機(jī)求解是指在5臺(tái)PC機(jī)上運(yùn)行分區(qū)并行算法，而單機(jī)求解是指在客戶機(jī)1上運(yùn)行該算法。從中我們可以明顯看出，使用區(qū)域分解法在MPI環(huán)境下并行求解燃燒器內(nèi)氣固兩相流動(dòng)問題有較高的加速比和并行效率。特別是其計(jì)算時(shí)間要低于同規(guī)模問題的單機(jī)串行計(jì)算時(shí)間約115 h，這在數(shù)值計(jì)算中是比較可觀的時(shí)間，可見在MPI并行平臺(tái)上對(duì)大型復(fù)雜問題進(jìn)行并行的數(shù)值求解有著很好的價(jià)值和前景。

表1 燃燒器內(nèi)分區(qū)并行計(jì)算并行效率

8 總結(jié)

本文在PC機(jī)群上通過MPI軟件包建立起了網(wǎng)絡(luò)并行計(jì)算平臺(tái)，并在此基礎(chǔ)上采用區(qū)域分解法對(duì)燃燒器一次風(fēng)噴咀和彎頭組件內(nèi)氣固兩相冷態(tài)流動(dòng)問題進(jìn)行了并行的數(shù)值求解。

本文的工作具體來說具有以下幾點(diǎn)意義：

(1)證實(shí)了在普通PC機(jī)集群上開展并行計(jì)算的可行性。

(2)本文自行開發(fā)出一套程序來模擬燃燒器一次風(fēng)噴咀和彎頭組件內(nèi)的氣固兩相流動(dòng)，正確預(yù)報(bào)其中兩相的速度、顆粒濃度及湍動(dòng)能分布。對(duì)于燃燒器結(jié)構(gòu)的改進(jìn)有一定的指導(dǎo)意義，也為以后進(jìn)行爐膛內(nèi)的數(shù)值模擬提供更加準(zhǔn)確可靠的邊界條件。

(3)通過區(qū)域分解方法實(shí)現(xiàn)分區(qū)并行求解，既減少了計(jì)算工作量，也提高了代數(shù)計(jì)算的品質(zhì)。

(4)在MPI并行平臺(tái)上成功實(shí)現(xiàn)對(duì)三維復(fù)雜流動(dòng)問題的并行求解，再一次證明了該并行算法具有良好、可靠和高效的工作性能。

(5)本文中所開發(fā)的基于MPI并行平臺(tái)的計(jì)算程序具有良好的數(shù)據(jù)通訊性和程序通用性，這為將來并行程序的模塊化打下了基礎(chǔ)。

總的說來，本論文的結(jié)果是符合流體力學(xué)的基本規(guī)律和實(shí)際情況的，使用區(qū)域分解的并行算法對(duì)模擬燃燒器一次風(fēng)噴咀和彎頭組件內(nèi)氣固兩相流問題也具有較高的效率。所以本文的工作對(duì)于在一般PC機(jī)集群上并行求解大規(guī)模復(fù)雜流動(dòng)問題是非常有參考和指導(dǎo)價(jià)值的。計(jì)算機(jī)數(shù)值模擬總是在朝著更高的計(jì)算精度、更好的區(qū)域適應(yīng)性和更強(qiáng)的健壯性方向發(fā)展。在計(jì)算機(jī)技術(shù)高速發(fā)展的今天，數(shù)值模擬對(duì)工程實(shí)際具有很重要的指導(dǎo)意義。

[1]馮俊凱，沈幼庭.鍋爐原理及計(jì)算 [M].第2版.北京:科學(xué)出版社，1992.

[2]呂濤，石濟(jì)民，林振寶.區(qū)域分解算法——偏微分方程數(shù)值解新技術(shù) [M].北京:科學(xué)出版社，1997.

[3]郁志輝.高性能計(jì)算并行編程技術(shù) [M].北京：清華大學(xué)出版社，2001.

[4]陸霄露.基于PC網(wǎng)絡(luò)求解復(fù)雜三維流動(dòng)的分區(qū)并行方法 [D].上海:上海交通大學(xué)，2003.

[5]陶文銓.計(jì)算傳熱學(xué)的近代進(jìn)展 [M].北京：科學(xué)出版社，2000.

[6]Patankar S V.郭寬良譯.傳熱和流體流動(dòng)的數(shù)值方法[M].北京:科學(xué)出版社,1984.

[7]蒲砢.關(guān)于氣固兩相流數(shù)值計(jì)算模型的探討 [J].重慶建筑高等?？茖W(xué)校學(xué)報(bào)，2000,10(4).

[8]周力行.湍流氣粒兩相流動(dòng)和燃燒的理論與數(shù)值模擬[M].北京：科學(xué)出版社，1994.

[9]陳國良.并行計(jì)算——結(jié)構(gòu)·算法·編程[M].北京:高等教育出版社，1999.

[10]陳國良.并行算法的設(shè)計(jì)和分析 [M].第2版.北京：高等教育出版社，2002.

[11]王小偉，郭力，唐德翔,等.兩相體系宏觀擬顆粒模擬并行算法 [J].化工學(xué)報(bào)，2004，55(5)：716-720.

[12]呂雄偉.基于網(wǎng)絡(luò)的循環(huán)流化床鍋爐整體數(shù)學(xué)模型的分布式算法研究 [D].重慶:重慶大學(xué)，2001.

[13]Brakkee E A,Segal A,Kassels C G M.A parallel domain decomposition algorithm for the incompressible Navier-Stokes equations[J].Simulation Practice and Theory,1995(3):363-389.

[14]呂曉斌，蘭黔章，朱自強(qiáng).求解Euler方程的區(qū)域分解方法與并行算法[J].計(jì)算物理，2000，17(4)：360-366.

[15]王開春，李樹民，朱國林.三維低速NS方程的并行計(jì)算 [J].空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20：88-93.