燃燒器內(nèi)氣固兩相流動的分區(qū)并行計算

符棟良

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

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

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

1 區(qū)域分解法

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

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

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

2.1 MPI介紹

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

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

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

2.2 MPI并行程序的主從模式

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

圖2 并行計算主從模式的流程

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

主從模式具有以下兩個優(yōu)點。

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

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

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

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

3 燃燒器的物理模型

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

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

圖3 帶有菱形體的一次風噴嘴及彎頭組件

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

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

在區(qū)域分解中，將整個計算區(qū)域分裂為六個子域Ω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)為其邊界，?Ω 為整個計算區(qū)域的邊界。圖4是燃燒器一次風組件整體的網(wǎng)格分布圖。其中在數(shù)值模擬中用到的爐膛部分的網(wǎng)格圖未給出，僅給出燃燒器一次風組件內(nèi)的網(wǎng)格圖。

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

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

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

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

6 并行方案

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

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

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

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

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

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

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

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

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

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

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

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

以下對Re=105時燃燒器一次風噴咀和彎頭組

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

7.1 流場分析

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

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

7.2 整流現(xiàn)象

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

7.3 湍動能分析

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

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

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

7.4 濃度場分析

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

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

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

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

7.5 并行性能分析

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

其計算公式如下：

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

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

tn——并行執(zhí)行時間，即程序在并行機群上的運行時間；

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

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

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

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

8 總結(jié)

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

本文的工作具體來說具有以下幾點意義：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[11]王小偉，郭力，唐德翔,等.兩相體系宏觀擬顆粒模擬并行算法 [J].化工學(xué)報，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]呂曉斌，蘭黔章，朱自強.求解Euler方程的區(qū)域分解方法與并行算法[J].計算物理，2000，17(4)：360-366.

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