基于最佳預(yù)后元模型的顆粒污垢特性全局敏感性分析

謝廣爍，張斯亮，何松，肖娟，王斯民

（1 西安交通大學(xué)化學(xué)工程與技術(shù)學(xué)院，陜西 西安 710049；2 上海藍(lán)濱石化設(shè)備有限責(zé)任公司，上海 201518）

換熱器廣泛應(yīng)用在化工、石油、電力、食品等工業(yè)生產(chǎn)過程中，是關(guān)系到國計民生的重要設(shè)備[1-4]。換熱管作為管殼式換熱器的主要元件，污垢積累會惡化傳熱性能，增加流動阻力，降低使用壽命。Akin? 等[5]發(fā)現(xiàn)煉油廠內(nèi)換熱器結(jié)垢造成的能源損失占總能耗的2%。顆粒污垢是換熱器內(nèi)主要的污垢類型之一，不僅會增加析晶污垢的成核數(shù)量，還會為生物的聚集、生長、繁殖提供場所，使多種污垢并存而難以去除[6-7]。

眾多學(xué)者分析了不同操作工況對顆粒沉積特性的影響[8-10]。劉洪濤等[11]研究了流動方式和流速對壁面上顆粒沉積速率的影響，發(fā)現(xiàn)水平流動下沉積速率隨顆粒量綱為1弛豫時間的增加而增加。張寧等[12]模擬了圓管內(nèi)的顆粒污垢生長特性，研究了流速、入口溫度和管壁粗糙度對污垢熱阻的影響。于曉燕[13]計算了圓管內(nèi)濃度、流速和溫度對顆粒沉積的影響，結(jié)果表明濃度和流速對污垢特性影響較為明顯。以上研究分析了常見操作參數(shù)對顆粒沉積的影響，如顆粒直徑、濃度、流速和溫度等，但并未對各參數(shù)的影響程度進(jìn)行定量描述。

敏感性分析作為定量評估模型或系統(tǒng)的輸入變量對輸出變量影響程度的重要方法，其中全局敏感性分析可分析每個變量及其交互作用對輸出變量的影響，量化設(shè)計變量的不確定性并有效識別關(guān)鍵參數(shù)，避免優(yōu)化無關(guān)參數(shù)導(dǎo)致有限的性能提升及計算負(fù)荷浪費[14-16]。Du 等[17]采用Taguchi 方法研究了螺旋角、管徑、管中心距等幾何參數(shù)對管殼式換熱器流動和傳熱性能的影響。Fesanghary 等[18]基于Sobol'全局敏感性量化了幾何參數(shù)對折流板換熱器總成本的影響，結(jié)果表明密封條數(shù)量是最不重要的參數(shù)，在優(yōu)化過程中將其剔除?；谠Ｐ偷拿舾行苑治鼍哂杏嬎懔啃　⑴c經(jīng)驗無關(guān)等優(yōu)點，得到廣泛應(yīng)用。元模型是一種表示輸入與輸出變量關(guān)系的近似數(shù)學(xué)算法，常用的元模型有多項式回歸、克里金和遺傳聚合模型等。由于工程問題中輸入與輸出變量關(guān)系的復(fù)雜性，元模型決定著敏感性分析的準(zhǔn)確性，因此最佳預(yù)后元模型（metamodel of optimal prognosis，MOP）應(yīng)運而生，從模型庫中自動選擇合適的元模型并處理變量，具有較高的通用性和準(zhǔn)確性[19]。Xiao 等[20]采用基于MOP 的全局敏感性分析，研究了折彎角、折彎率、相對高度和殼程流速對旋梯式螺旋折流板換熱器流動傳熱性能的影響，結(jié)果表明殼側(cè)流速和折彎角對單位壓降傳熱系數(shù)影響較為明顯，總影響分別為70.34%、25.01%。

內(nèi)插件不僅增強管內(nèi)傳熱，還具備在線防污除垢的作用，根據(jù)不同的運動方式分為固定式和轉(zhuǎn)動式。自轉(zhuǎn)式內(nèi)插件由于流體流過產(chǎn)生動力矩使其旋轉(zhuǎn)，管內(nèi)壓降通常低于靜止內(nèi)插件的壓降[21]，且比電機帶動下的主動式旋轉(zhuǎn)內(nèi)插件減少了額外動力消耗。因此，本文以內(nèi)置轉(zhuǎn)子防結(jié)垢為研究背景，基于MOP 對傳熱管內(nèi)顆粒污垢特性進(jìn)行全局敏感性分析。顆粒直徑、入口流速、顆粒濃度和入口溫度為輸入變量，沉積率和污垢熱阻為輸出參數(shù)，并且分析對擴散沉積、湍流泳沉積、熱泳沉積和重力沉積四個主要沉積速率的影響。訓(xùn)練數(shù)據(jù)空間通過數(shù)值計算獲得，采用了歐拉-歐拉模型和顆粒污垢模型，得到了輸入?yún)?shù)對輸出參數(shù)的主影響和交互影響，為傳熱管內(nèi)防結(jié)垢找尋關(guān)鍵變量。

1 數(shù)值計算方法與幾何模型

1.1 幾何模型

如圖1所示，傳熱管內(nèi)徑28mm，長度300mm；轉(zhuǎn)子由一個六邊形面扭轉(zhuǎn)90°得到，并在徑向最大處兩端各倒出一個半徑10mm 的圓角。轉(zhuǎn)子直徑26mm，節(jié)距78mm，為了方便轉(zhuǎn)子在實際應(yīng)用中的安裝，其左右兩端各加工一個直徑5mm、高10mm的圓臺。

圖1 幾何模型(單位：mm)

1.2 數(shù)學(xué)模型

1.2.1 控制方程歐拉-歐拉模型通過壓力和相間交換系數(shù)實現(xiàn)動量方程與連續(xù)性方程的耦合。對于泥漿流動、水力輸送和顆粒沉積的計算，使用歐拉-歐拉模型較為合適[22]。計算中進(jìn)行了以下簡化假設(shè)：①流動和傳熱都是穩(wěn)態(tài)過程；②忽略換熱管向環(huán)境的熱量損失；③不考慮質(zhì)量力和流體的黏性耗散。因此，歐拉-歐拉模型的控制方程如下式(1)～式(3)。

連續(xù)性方程

本文采用RNGk-ε模型進(jìn)行湍流求解，其加入了旋流對湍流的影響，提高了預(yù)測旋流和快速應(yīng)變流的精度[23]，如式(4)、式(5)。

湍動能k方程

1.2.2 顆粒污垢模型

顆粒污垢特性是沉積和剝蝕迭加后的結(jié)果，學(xué)者們采用不同的數(shù)值模型研究了顆粒沉積過程，如拉格朗日粒子跟蹤法、歐拉兩相法、離散粒子模型（DPM）和Mixture模型等[24-27]。文獻(xiàn)[28]詳細(xì)介紹了本文采用的顆粒污垢模型，考慮了水中微米顆粒沉積的四種有效沉積機制：擴散沉積、湍流泳沉積、熱泳沉積和重力沉積。顆粒沉積率計算如式(6)和式(7)。

顆粒剝蝕受壁面剪切力和污垢黏結(jié)強度因子影響，剝蝕率計算如式(8)。

假定污垢均勻分布在換熱管內(nèi)壁面上，且污垢成分一致，則污垢熱阻如式(9)。

1.3 邊界條件

重力加速度為y軸負(fù)方向，連續(xù)相為水，分散相為氧化鎂顆粒，氧化鎂顆粒的直徑是20～50μm，顆粒濃度為0.4～4.0kg/m3。轉(zhuǎn)子和管壁的材質(zhì)分別為鋁和鐵，水的黏度隨溫度變化較大（見表1），其余物性參數(shù)均為常物性（見表2）。采用速度入口和壓力出口邊界，入口速度為0.15～0.35m/s，介質(zhì)入口溫度為293～323K，管壁溫度恒定為343K。采用多重參考系法模擬轉(zhuǎn)子的轉(zhuǎn)動，將流體域劃分為轉(zhuǎn)子所在的旋轉(zhuǎn)域及其他靜止域。

表1 不同溫度下水的黏度[29]

表2 物性參數(shù)

1.4 全局敏感性分析方法

假設(shè)輸入與輸出之間的關(guān)系可以描述為Y=f(x)，其中Y=(y1,…,ym)是標(biāo)量輸出函數(shù)，x=(x1, … ,xn)是n維自由輸入變量。根據(jù)Sobol'[30]提出的全局敏感性指標(biāo)，一階效應(yīng)敏感性，也稱為主影響，用來衡量xi對y的單獨貢獻(xiàn)，如式(10)。

式中，V(Y)為f(x)的無條件方差；V[E(Y|xi)]為fi(xi)的方差。

由于一階敏感性僅表現(xiàn)為xi的變化，只有各個輸入變量的解耦效應(yīng)，而高階敏感性指標(biāo)通常不計算。Homma 和Saltelli[31]提出將x=(x1,…,xn)分成兩部分，即xi和其互補子集x～i，互補子集包括除xi以外的所有因素，總影響如式(11)。

STi=Si+Si,～i= 1 -S～i(11)

式中，Si,～i為xi與其他所有元素的相互影響；S～i為x～i的主影響且不包括xi對應(yīng)的任何影響。

最佳預(yù)后元模型的建立是一個自動選擇過程，元模型庫包括多項式回歸、最小二乘法和克里金模型。此外，考慮Box-Cox變換以提高元模型的近似質(zhì)量，該冪變換方法能夠穩(wěn)定方差，增強正態(tài)分布，減少非線性效應(yīng)。變量篩選一方面是采用由皮爾遜積矩相關(guān)系數(shù)決定的顯著性濾波器，只有當(dāng)輸入變量之間的相關(guān)誤差在95%～99%以內(nèi)時才考慮其對輸出變量的顯著性。另一方面采用重要性濾波器，其定義為全模型與約簡模型的調(diào)整決定系數(shù)（adj.R2）之差。CoImin的最小閾值在1%～5%之間，如果CoI不小于CoImin，則將考慮該變量。在交叉驗證的基礎(chǔ)上， 提出預(yù)后系數(shù)（coefficient of prognosis，CoP）衡量MOP 的近似質(zhì)量。結(jié)合基于方差的敏感性指標(biāo)，將CoP 與MOP 確定的總影響敏感性指標(biāo)相乘，計算出單個輸入變量的方差貢獻(xiàn)，如式(12)～式(15)。

圖2 展示了基于MOP 的全局敏感性分析流程。輸入變量為顆粒直徑d（20～50μm）、顆粒濃 度c（0.4～4.0kg/m3）、入 口 速 度v（0.15～0.35m/s）、介質(zhì)入口溫度T（293～323K），沉積率和污垢熱阻作為輸出參數(shù)，采用增強的旋轉(zhuǎn)中心差分方法生成49 個實驗點。

圖2 基于MOP的全局敏感性分析

2 網(wǎng)格無關(guān)性驗證與模型驗證

當(dāng)顆粒直徑為20μm、顆粒濃度為0.4kg/m3、入口溫度為303K、入口流速為0.3m/s 時，對傳熱管模型進(jìn)行網(wǎng)格獨立性驗證。流體域采用四面體網(wǎng)格，并局部加密管壁邊界層和轉(zhuǎn)子壁面網(wǎng)格。通過改變整體網(wǎng)格、管壁第一層網(wǎng)格和轉(zhuǎn)子壁面網(wǎng)格尺寸來調(diào)整網(wǎng)格數(shù)量，結(jié)果如圖3所示。隨著網(wǎng)格數(shù)量的增加，傳熱管內(nèi)的沉積率和壓降都是先減小后趨于穩(wěn)定。當(dāng)網(wǎng)格數(shù)量為1727959 時，沉積率為2.62×10-4kg/(m2·s)，相較于網(wǎng)格數(shù)為1250704 時變化了2.2%，同時壓降變化了1.4%。綜合考慮計算精度和計算時間，本文采用的網(wǎng)格數(shù)為1250704。

圖3 網(wǎng)格獨立性測試

通過復(fù)現(xiàn)文獻(xiàn)[32]的實驗結(jié)果，對MRF模型和RNGk-ε模型計算內(nèi)置轉(zhuǎn)子傳熱管內(nèi)流場和溫度場的準(zhǔn)確性進(jìn)行驗證。以內(nèi)管中的流體為研究對象，入口溫度為313K，管壁溫度為295K，傳熱管內(nèi)徑為20mm， 扭帶寬度和長度分別為14mm、2000mm，節(jié)距比為6.0，建立了與文獻(xiàn)中一致的物理模型，模擬結(jié)果與實驗數(shù)據(jù)對比如圖4所示，模擬值與實驗值的變化趨勢一致，努塞爾數(shù)的偏差最小為6.1%，平均偏差為8.4%，表明內(nèi)置轉(zhuǎn)子傳熱管內(nèi)流場和溫度場計算的精度較高。

圖4 內(nèi)置轉(zhuǎn)子傳熱管Nu數(shù)計算驗證

通過復(fù)現(xiàn)文獻(xiàn)[33]中不同入口流速下圓管內(nèi)顆粒污垢熱阻測試實驗，對本文所采用的顆粒污垢數(shù)學(xué)模型進(jìn)行驗證，計算結(jié)果與實驗結(jié)果對比如圖5 所示?？梢园l(fā)現(xiàn)，隨著時間的增加污垢熱阻逐漸增加，并在一段時間后達(dá)到漸近值。模擬值與實驗值增長速度相近，三種工況下漸近值的相對偏差均在12.0%以內(nèi)，最小偏差為6.2%，符合兩相流計算要求。漸近值最能反映換熱管長期運行時所產(chǎn)生的污垢熱阻，因此本文后面提及的污垢熱阻均是漸近值。

圖5 污垢熱阻模擬值與實驗數(shù)據(jù)對比曲線

圖6 四個沉積速率的總影響CoP矩陣

3 結(jié)果與討論

3.1 沉積速率全局敏感性分析

通過構(gòu)建MOP 實現(xiàn)最重要參數(shù)和最佳元模型的自動檢測，并驗證了預(yù)后質(zhì)量。由于對沉積速率進(jìn)行量綱為1化時消除了顆粒濃度的影響，所以本節(jié)主要研究入口溫度、入口流速和顆粒直徑三個因素對各量綱為1沉積速率的影響。

3.2 沉積率與污垢熱阻全局敏感性分析

圖7 示出了四個輸入變量對沉積率和污垢熱阻影響的CoP 矩陣。從圖中可以看出，md和Rf的全模型CoP 分別為97.8%、96.2%，預(yù)后模型都是基于BoX-CoX 變換的一階多項式回歸。單個CoP的最大值均為d所對應(yīng)數(shù)值，說明d是影響顆粒沉積特性的最顯著因素。對于md，d和v的總影響較為明顯，c和T的總影響接近且較小。對于Rf，除d外其余三個變量的總影響相近，且遠(yuǎn)小于d的總影響。

圖7 沉積率和污垢熱阻的總影響CoP矩陣

圖8 和圖9 示出了四個變量對md、Rf的主影響和交互影響柱狀圖。圖8 可以看出，四個變量對md的主影響和總影響的排列順序一致，d對md的主影響比交互影響高出了71.6%，而v、c和T對md的主影響和交互影響十分接近。主影響和交互影響的總和分別為73.0%、59.7%，說明各變量之間的交互作用對md有著重要的影響。圖9 可以發(fā)現(xiàn)，d對Rf的主影響明顯高于交互影響，主影響為47.3%，而交互影響為12.8%，v和c對Rf的主影響均略高于交互影響，T對Rf的主影響比交互影響低48.1%。主影響和交互影響的總和分別為75.7%、42.5%，各變量之間的交互作用依然不可忽視。

圖8 md的主影響和交互影響柱狀圖

圖9 Rf的主影響和交互影響柱狀圖

圖10和圖11分別示出了混合項對md和Rf的預(yù)測響應(yīng)，每個云圖均在其他參數(shù)取變化范圍中間值時得到。由圖10 可知，最大的沉積率都出現(xiàn)在云圖的右上角，說明隨著四個輸入變量的增加，沉積率增大。隨著c的提升，更多的顆粒與壁面發(fā)生碰撞，單位時間內(nèi)沉積到壁面上的顆粒增加。熱泳運動中顆粒由溫度較高向溫度較低的地方移動，隨著T的增加，壁面與流體的溫差減小，熱泳運動對顆粒沉積的抑制作用減弱；同時T提高使流體的黏度減小，湍流泳沉積迅速增加。d越大，重力的作用越明顯，重力沉積快速增加。v越大，單位時間內(nèi)流體攜帶更多的顆粒流向壁面，沉積率增加。

圖10 預(yù)測混合項對md的響應(yīng)

圖10(a)示出了c和T在v=0.25m/s、d=35μm 時對md的影響，可以看出，在參數(shù)變化范圍內(nèi)沉積率一直偏低，且沉積率等值線與豎直方向呈45°夾角，說明c和T對md的綜合影響較小，且兩者對md的影響效果相近。在圖10(b)、(c)、(f)中，當(dāng)d=20μm 時，c、T、v對md的影響很?。浑S著d的增加，c、T、v對md的影響逐漸增加，說明c、T、v對md的局部敏感度逐漸增加。同理，在圖10(a)、(d)、(e)中，等值線從圖中左下角到右上角越來越密集，說明隨著各變量的增加，輸入變量對md的局部敏感度越來越大。

圖11 表明污垢熱阻的最小值都出現(xiàn)在每個云圖的左下角，隨著c、T、d、v的增加污垢熱阻逐漸增大，與沉積率的變化規(guī)律一致，且圖11(a)與圖10(a)類似，表明c和T對Rf的影響效果相近。在圖11(b)、(c)、(f)中出現(xiàn)了較高污垢熱阻（紅色區(qū)域），說明d和c、d和T、d和v對Rf的綜合影響十分顯著。保持c=2.2kg/m3、T=308K 不變，當(dāng)d=20～26μm，Rf幾乎不隨v變化。此外，等值線從左下角到右上角逐漸變得密集，表示輸入變量對Rf的影響越來越明顯，表明隨輸入變量的增加，各輸入變量對Rf的局部敏感度越來越高。

4 結(jié)論

為了探究影響內(nèi)置轉(zhuǎn)子傳熱管內(nèi)顆粒污垢的主要因素，采用基于MOP 的全局敏感性分析方法，研究了顆粒直徑、入口流速、顆粒濃度和入口溫度四個因素對顆粒污垢特性的影響，并區(qū)分了不同因素的影響程度。

（2）沉積率和污垢熱阻的最佳預(yù)后元模型均為基于BoX-CoX 變換的一階多項式回歸，相應(yīng)的全模型CoP分別為97.8%、96.2%。

（3）顆粒直徑是對沉積率和污垢熱阻影響最顯著的因素，總影響分別為52.7%、60.2%。入口流速對沉積率也有較明顯的影響，總影響為41.9%；但入口流速、顆粒濃度和入口溫度對污垢熱阻的影響明顯小于顆粒直徑。

（4）變量之間的交互作用對沉積率和污垢熱阻的影響不可忽視，四個變量對沉積率交互影響的總和為59.7%，對污垢熱阻交互影響的總和為42.5%，入口溫度對污垢熱阻的主影響比交互影響低48.1%。

符號說明

c—— 顆粒濃度，kg/m3

d—— 顆粒粒徑，μm

H—— 焓值，J

l—— 管子長度，mm

md—— 沉積率，kg/(m2·s)

mf—— 顆粒凈沉積總量，kg/m2

mr—— 剝蝕率，kg/(m2·s)

p—— 壓強，Pa

Rf—— 污垢熱阻，m2·K/W

T—— 入口溫度，K

u—— 量綱為1沉積速率

u—— 量綱為1重力沉積速率

u—— 量綱為1擴散沉積速率

u—— 量綱為1湍流泳沉積速率

u—— 量綱為1熱泳沉積速率

u*—— 壁面剪切速度，m/s

v—— 入口流速，m/s

α—— 體積分?jǐn)?shù)

λf—— 污垢層熱導(dǎo)率，W/(m·K)

μ—— 動力黏度，Pa·s

ξ—— 污垢黏結(jié)強度因子，N·s/m2

ρ—— 密度，kg/m3

τw—— 剪切力，N/m2

下角標(biāo)

d—— 沉積

f—— 污垢

i,j—— 矢量分量

k—— 狀態(tài)相，k=l時代表流體相，k=p時代表顆粒相