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

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

（1 西安交通大學化學工程與技術學院，陜西 西安 710049；2 上海藍濱石化設備有限責任公司，上海 201518）

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

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

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

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

1 數值計算方法與幾何模型

1.1 幾何模型

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

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

1.2 數學模型

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

連續性方程

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

湍動能k方程

1.2.2 顆粒污垢模型

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

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

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

1.3 邊界條件

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

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

表2 物性參數

1.4 全局敏感性分析方法

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

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

由于一階敏感性僅表現為xi的變化，只有各個輸入變量的解耦效應，而高階敏感性指標通常不計算。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對應的任何影響。

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

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

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

2 網格無關性驗證與模型驗證

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

圖3 網格獨立性測試

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

圖4 內置轉子傳熱管Nu數計算驗證

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

圖5 污垢熱阻模擬值與實驗數據對比曲線

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

3 結果與討論

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

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

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

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

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

圖8 和圖9 示出了四個變量對md、Rf的主影響和交互影響柱狀圖。圖8 可以看出，四個變量對md的主影響和總影響的排列順序一致，d對md的主影響比交互影響高出了71.6%，而v、c和T對md的主影響和交互影響十分接近。主影響和交互影響的總和分別為73.0%、59.7%，說明各變量之間的交互作用對md有著重要的影響。圖9 可以發現，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的預測響應，每個云圖均在其他參數取變化范圍中間值時得到。由圖10 可知，最大的沉積率都出現在云圖的右上角，說明隨著四個輸入變量的增加，沉積率增大。隨著c的提升，更多的顆粒與壁面發生碰撞，單位時間內沉積到壁面上的顆粒增加。熱泳運動中顆粒由溫度較高向溫度較低的地方移動，隨著T的增加，壁面與流體的溫差減小，熱泳運動對顆粒沉積的抑制作用減弱；同時T提高使流體的黏度減小，湍流泳沉積迅速增加。d越大，重力的作用越明顯，重力沉積快速增加。v越大，單位時間內流體攜帶更多的顆粒流向壁面，沉積率增加。

圖10 預測混合項對md的響應

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

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

4 結論

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

（2）沉積率和污垢熱阻的最佳預后元模型均為基于BoX-CoX 變換的一階多項式回歸，相應的全模型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

α—— 體積分數

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

μ—— 動力黏度，Pa·s

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

ρ—— 密度，kg/m3

τw—— 剪切力，N/m2

下角標

d—— 沉積

f—— 污垢

i,j—— 矢量分量

k—— 狀態相，k=l時代表流體相，k=p時代表顆粒相