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

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

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

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

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

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

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

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

1.1 幾何模型

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

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

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

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

連續(xù)性方程

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

湍動(dòng)能k方程

1.2.2 顆粒污垢模型

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

顆粒剝蝕受壁面剪切力和污垢黏結(jié)強(qiáng)度因子影響，剝蝕率計(jì)算如式(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)動(dòng)，將流體域劃分為轉(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對(duì)y的單獨(dú)貢獻(xiàn)，如式(10)。

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

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

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

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

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

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

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

當(dāng)顆粒直徑為20μm、顆粒濃度為0.4kg/m3、入口溫度為303K、入口流速為0.3m/s 時(shí)，對(duì)傳熱管模型進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證。流體域采用四面體網(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 時(shí)，沉積率為2.62×10-4kg/(m2·s)，相較于網(wǎng)格數(shù)為1250704 時(shí)變化了2.2%，同時(shí)壓降變化了1.4%。綜合考慮計(jì)算精度和計(jì)算時(shí)間，本文采用的網(wǎng)格數(shù)為1250704。

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

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

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

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

圖5 污垢熱阻模擬值與實(shí)驗(yàn)數(shù)據(jù)對(duì)比曲線

圖6 四個(gè)沉積速率的總影響CoP矩陣

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

符號(hào)說明

c—— 顆粒濃度，kg/m3

d—— 顆粒粒徑，μm

H—— 焓值，J

l—— 管子長度，mm

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

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

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

p—— 壓強(qiáng)，Pa

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

T—— 入口溫度，K

u—— 量綱為1沉積速率

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

u—— 量綱為1擴(kuò)散沉積速率

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

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

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

v—— 入口流速，m/s

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

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

μ—— 動(dòng)力黏度，Pa·s

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

ρ—— 密度，kg/m3

τw—— 剪切力，N/m2

下角標(biāo)

d—— 沉積

f—— 污垢

i,j—— 矢量分量

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