王才謙 劉石

1.華北電力大學(xué)能源動力與機(jī)械工程學(xué)院；2.華北電力大學(xué)控制與計算機(jī)工程學(xué)院

針對流場中存在復(fù)雜機(jī)械結(jié)構(gòu)而導(dǎo)致得濃度難以測量的問題，提出了基于塔克分解的濃度重建算法。采用塔克分解對樣本數(shù)據(jù)庫進(jìn)行特征提取，將其運用于濃度分布重建。為了驗證算法的可行性，進(jìn)行了泄漏仿真實驗。仿真實驗結(jié)果表明濃度重建算法可以快速重建濃度分布。最后將本方法與基于主成分分析的重建方法進(jìn)行比較，實驗數(shù)據(jù)表明本方法在三維濃度重建上具有更好的重建精度。

流場中的濃度信息可以通過實驗測量和數(shù)值計算獲得。其中實驗測量主要可分為接觸式與非接觸式兩大類。接觸式測量通過將傳感器放置在流場中進(jìn)行采樣測量，只能獲取流場少數(shù)位置的濃度值，無法獲取全場的瞬時濃度。若要進(jìn)行多點同時采樣，需要大量的采樣儀器與分析設(shè)備，測量成本高。非接觸式測量利用流體的光學(xué)[1]或者聲學(xué)[2]特性進(jìn)行濃度測量，但是這些設(shè)備也會有結(jié)構(gòu)復(fù)雜、布置困難、易受環(huán)境影響等問題，一定程度上制約了非接觸時測量的應(yīng)用。近年來，隨著計算流體力學(xué)的發(fā)展與計算機(jī)水平的進(jìn)步，數(shù)值計算廣泛地運用于流場分析中。但是要想求解流場的濃度值信息，需要獲取流場的邊界條件，因此在邊界條件未知時，就難以使用數(shù)值計算方法。而且使用數(shù)值計算方法需要大量的計算資源和時間，很難實時獲取流場的濃度值信息。

本文的主要工作是將數(shù)值計算方法和測量方法結(jié)合起來，提出一種基于塔克分解的濃度重建算法。該算法將CFD 數(shù)值模擬得到的流場濃度值作為先驗信息，然后將先驗信息降維提取主要特征后，結(jié)合測量的少量濃度值，重建出全場的濃度分布。

1 算法原理

基于特征提取的數(shù)據(jù)處理方法最早運用于屋蓋風(fēng)壓場的重建[3]，后來又廣泛運用于其他流場參數(shù)的重建過程中：QIN 等人[4]將主成分分析、CFD 計算數(shù)據(jù)結(jié)合起來，實現(xiàn)了對風(fēng)場分布的快速預(yù)測；Sun 等人[5]提出一種基于本征正交分解的風(fēng)速重建方法，并且對測量點位置進(jìn)行了優(yōu)化，減少了重建誤差；陳敏鑫等人[6]將主成分分析引入到溫度分布重建工作中，利用經(jīng)驗數(shù)據(jù)和少量測量數(shù)據(jù)快速準(zhǔn)確地完成了溫度分布重建。但是以上方法研究的均是二維平面上流場參數(shù)的重建，所使用的均是基于向量或矩陣的數(shù)據(jù)降維方法。如果將這些方法運用于三維流場的重建，會在一定程度上造成三維空間特征的損失，為了避免這一情況出現(xiàn)，本文引入了塔克分解作為降維算法對先驗信息進(jìn)行特征提取。

1.1 塔克分解

隨著對數(shù)據(jù)的深入研究，矩陣已經(jīng)很難完整表述出數(shù)據(jù)的時空關(guān)系。而張量作為可以表征高階數(shù)據(jù)的多維數(shù)組，可以完整地保留各維度的信息。張量的構(gòu)成如圖1所示。

圖1 張量示意圖Fig.1 Tensor diagram

1963 年Tucker 首次提出了塔克分解，將一個n 階高維的張量，得到n 個分解因子矩陣和一個n 階低維的核心張量。n 階高維張量的分解過程的公式如式（1）所示：

圖2 三階張量塔克分解示意圖Fig.2 Schematic diagram of the third-order tensor Tucker decomposition

1.2 濃度分布重建算法

利用經(jīng)驗數(shù)據(jù)或者測量數(shù)據(jù)獲得n 個工況下的三維濃度分布，組成一份先驗數(shù)據(jù)集XT={X1，X2，…，。對先驗數(shù)據(jù)集XT 進(jìn)行塔克分解將其分解成核心張量與分解因子矩陣模態(tài)積的形式，如式（2）所示。

將XT 與F(4)內(nèi)的元素展開，公式（2）可轉(zhuǎn)化為如式（3）所示形式：

式中，bi為4 模態(tài)方向上的分解因子矩陣中的第i 個因子向量，

由公式（3）可得：

那么對于待重建的三維濃度分布XR,如式（5）所示：

式中，b 為4 模態(tài)方向上待重建的三維濃度分布XR所對應(yīng)的分解因子向量。

根據(jù)模態(tài)積的定義將式（5）從張量形式轉(zhuǎn)化為矩陣與向量相乘的形式，如式（6）所示：

式中，xR為將三維濃度分布XR中的元素轉(zhuǎn)換為向量形式，維度為(n1×n2×n3)×1；Aunfold為將重建系數(shù)張量沿著模態(tài)4 方向展開的矩陣形式，維度為(n1×n2×n3)×m4。

少量的測量數(shù)據(jù)與整個流場的濃度分布的位置關(guān)系可由式（7）表示。

式中，y 為實際測量得到的數(shù)據(jù)，維數(shù)為r×1，其中r為測點數(shù)目；M 為濃度測量矩陣，維數(shù)為r×(n1×n2×n3)。濃度測量矩陣M 為0-1 矩陣，每一行只有一個元素的值為1，其余元素的值為0。在第n 行（n 取1 到r），值為1 的元素所在的位置，代表了第n 個測量點與全部濃度點的相對位置。

為了使重建的三維濃度分布與實際的三維濃度分布相等，即x=xR，把式（6）代入到式（7）中，可以得到如式（8）所示：

在式（8）中，y 通過測量獲得，M 為測量矩陣，在決定好測點位置后為已知量，通過對樣本數(shù)據(jù)張量XT 進(jìn)行塔克分解后獲得，因此可以先用最小二乘法求解出式（8）中的b，再代入式（6）即可求解出待重建的三維濃度分布。

最后為了檢驗重建效果，定義相對重建誤差：

基于塔克分解的濃度重建算法如圖3 所示。從圖3中可以看出，完成重建算法的第一步是利用CFD 進(jìn)行數(shù)值模擬，得到邊界條件下濃度場分布數(shù)據(jù)集；第二步是對樣本張量進(jìn)行塔克分解，求解其核心張量與各個模態(tài)方向上的分解因子矩陣；第三步是求解重建系數(shù)張量，并將其展開為重建系數(shù)矩陣；第四步是得到各測點的位置與測量的濃度值；最后計算待重建的濃度分布向量所對應(yīng)的分解因子向量b，進(jìn)而求出重建結(jié)果。

圖3 基于塔克分解的濃度重建算法流程圖Fig.3 Flow chart of concentration reconstruction algorithm based on tak decomposition

2 案例驗證

為了驗證算法的可靠性，建立泄漏模型。使用前處理軟件ICEM 建立幾何結(jié)構(gòu)并劃分網(wǎng)格，使用Fluent17.0進(jìn)行仿真求解。

2.1 模型幾何信息

泄漏模型如圖4所示，整個容器的計算尺寸為200mm×200mm×875mm，水流入口取水在容器中經(jīng)過穩(wěn)流裝置后的截面，金屬棒下端固定在距水流入口50mm 高的位置上，每個金屬棒長565mm，為了簡化計算，固定裝置忽略不計，水流出口的截面尺寸為60mm×40mm。冷卻水從下方的水流入口進(jìn)入，經(jīng)過金屬棒然后從水流出口流出。中間的金屬棒在高度為420mm 的位置上存在一個直徑為3mm 的圓形破口，方向朝著x 軸正方向。泄漏物會由于濃度差以及流體攜帶等機(jī)制向水中擴(kuò)散和傳輸，最終隨水流從出口排出。

圖4 流動模型結(jié)構(gòu)圖Fig.4 Flow model structure diagram

從圖4 中可以看出，9 根金屬棒遮擋了整個容器下半部分的濃度分布，很難直接測量這些區(qū)域的濃度值，而金屬棒上方的濃度分布是可以利用前文中提到的測量方式測量的。因此泄漏模型驗證的目的是利用金屬棒上方區(qū)域可測量的濃度分布來重建出金屬棒遮擋區(qū)域不可測量的濃度分布。

2.2 網(wǎng)格劃分

結(jié)構(gòu)化網(wǎng)格具有生成速度快、網(wǎng)格質(zhì)量高、數(shù)據(jù)結(jié)構(gòu)簡單等優(yōu)點，但是對于不規(guī)則區(qū)域的適應(yīng)性較差。由于前文中對本模型進(jìn)行了簡化，幾何結(jié)構(gòu)并不復(fù)雜，因此使用結(jié)構(gòu)化網(wǎng)格進(jìn)行網(wǎng)格劃分。

一般來說網(wǎng)格數(shù)量越多，計算結(jié)果就能更加精確，但是計算耗時也會變長。為了在從稀疏到密集的5 組網(wǎng)格中選取一個合適的網(wǎng)格節(jié)點數(shù)，在保證精度的前提下減少計算耗時，對這5 組網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性檢驗。5組網(wǎng)格的網(wǎng)格數(shù)分別為1131852、1729692、2658972、4526216、6142286。選取中間金屬棒頂端中心到水流出口中心所連成的線段上100 個點的流速作為檢驗計算結(jié)果的判斷標(biāo)準(zhǔn)，不同網(wǎng)格數(shù)目下樣本點流速如圖5 所示。

圖5 不同網(wǎng)格數(shù)下各樣本點流速Fig.5 Flow velocity of each sample point under different grid numbers

觀察樣本點折線可得，使用以上網(wǎng)格所計算出的結(jié)果差別不大，考慮到在網(wǎng)格數(shù)為1131852 時計算速度最快，因此在本文中仿真使用網(wǎng)格數(shù)為1131852 的網(wǎng)格。

2.3 計算模型

泄漏實驗的泄漏物為羅丹明B 水溶液，由于溶液的濃度不高，且不會與水發(fā)生化學(xué)反應(yīng)，其物理性質(zhì)仍按照水的性質(zhì)設(shè)置。泄漏物與水都視為不可壓縮流體，對于水流的湍流流動，采用k-ε 湍流模型方程并結(jié)合連續(xù)性方程和動量守恒方程進(jìn)行三維仿真。選擇速度、壓力的耦合求解算法，動量和湍動能采用二階迎風(fēng)格式求解。水流入口和泄漏口設(shè)定為速度入口，水流出口設(shè)定為壓力出口，出口壓力設(shè)為0pa。殘差的收斂條件設(shè)為10e-4，根據(jù)殘差和出口、入口質(zhì)量守恒來判斷計算是否收斂。

3 重建結(jié)果

本文將重建區(qū)域的三維空間用幾個堆疊的平面來表示。由于計算結(jié)果中越靠近z=0m 處的濃度值較高，越靠近z=±0.1m 處的濃度值較低，而濃度值較低的點對濃度重建的意義不大，因此本文中選取z=-0.04m，z=-0.02m，z=0m，z=0.02m，z=0.04m 五個平面來代表整個三維流體域，如圖6 陰影部分所示。

圖6 堆疊平面示意圖Fig.6 Stacked floor plan

每個平面中濃度點的排列方式為：在x 方向上從x=-0.1m 開始到x=0.1m 結(jié)束，每隔0.002m 設(shè)置一個濃度點，共設(shè)置101 個；在y 方向上從y=0.42m 開始到y(tǒng)=0.82m 結(jié)束，每隔0.002m 設(shè)置一個濃度點，共設(shè)置201個，每一個平面的濃度點數(shù)目為20301 個。本文根據(jù)不同的水流速與泄漏物流速設(shè)置36 個邊界條件作為樣本工況，水流速分別為0.01、0.02、0.03、0.04、0.05 和0.06m/s，泄漏物流速分別為0.01、0.02、0.03、0.04、0.05 和0.06m/s。將樣本工況計算完成后，建立一個101×201×5×36 的四維張量樣本數(shù)據(jù)庫。

另外，在實際測量的過程中，只知道測量點的濃度分布，不知道水流速度與泄漏物泄漏速度這些邊界條件，測量的工況基本上都是在樣本工況之外。因此，想要驗證濃度重建算法的實用性，還要驗證待重建的濃度分布工況位于樣本工況外的重建結(jié)果，因此除了36 個樣本工況外還設(shè)置了5 個測試工況，構(gòu)成了一個維數(shù)為101×201×5×5 的四維張量測試數(shù)據(jù)集。

從高度在0.62m 以上的濃度點（即沒有被金屬棒遮擋的區(qū)域的濃度點）中隨機(jī)選取400 個點作為測量點數(shù)據(jù)帶入到重建算法中，三維濃度重建所使用的測試工況如表1 所示。

表1 測試工況參數(shù)Tab.1 Test condition parameters

經(jīng)過濃度重建后，5 個測試工況的重建誤差如表2所示。

表2 測試工況重建誤差Tab.2 Reconstruction error of test conditions

從表2 可得，各測試工況的重建誤差均在10%以下，重建時間在1s內(nèi)，因此只要能夠及時獲取測量點的濃度值，就可以獲得全場的瞬時濃度。重建所得測試工況5的濃度分布如圖7 所示，由于各平面的濃度值的數(shù)量級差距較大，故將各平面的濃度分布分開展示。

由圖7 可以看出，在測試工況不在樣本工況范圍內(nèi)的情況下，使用基于塔克分解的濃度重建分布算法在z=±0.02m 和z=0m 這三個平面有較好的重建結(jié)果。對于z=±0.04m 兩個平面來說，雖然在濃度分布圖上的數(shù)值計算和重建的結(jié)果差距較大，但是這兩個平面的濃度值很低，為10-9的數(shù)量級，雖然相對重建誤差較大，但是絕對誤差仍然在一個較小的水平。因此可以認(rèn)為本算法可以實現(xiàn)對模型濃度分布的有效重建。

圖7 各平面數(shù)值計算與重建濃度分布對比Fig.7 Comparison of numerical calculation and reconstructed concentration distribution in each plane

此外，在三維濃度分布重建算法的研究中，可以將一個工況下的三維濃度點進(jìn)行編號，將其儲存在一個向量中，計算多個工況的濃度值就形成一個樣本矩陣，再采用基于主成分分析的方法來提取特征向量[7]。為了比較該方法與基于塔克分解方法的優(yōu)劣，另外采用了基于主成分分析的濃度重建算法與基于塔克分解的重建算法進(jìn)行對比。

如圖8 所示，相較于基于塔克分解的重建算法，使用基于主成分分析的重建算法的重建誤差在各測試工況下都要更高。而且使用主成分分析的過程中對先驗信息的處理會涉及到協(xié)方差矩陣的計算，在先驗信息維數(shù)較高時，計算機(jī)甚至?xí)霈F(xiàn)內(nèi)存不足無法處理的情況。因此使用塔克分解的重建算法相較于使用主成分分析的重建算法在三維濃度重建方面具有更好的效果。

圖8 各工況下基于主成分分析與塔克分解的重建算法的重建誤差對比Fig.8 Comparison of reconstruction errors of reconstruction algorithms based on principal component analysis and Tucker decomposition under various working conditions

4 結(jié)語

本文針對由于技術(shù)條件和環(huán)境影響難以獲取流體域內(nèi)部濃度信息的問題，將特征提取引入濃度分布重建的研究中，結(jié)合了實驗測量與數(shù)值計算方法，提出了基于塔克分解的三維濃度重建算法，并從原理上論證了算法的可行性，并總結(jié)了進(jìn)行濃度分布重建的計算流程。然后對于泄漏過程進(jìn)行數(shù)值模擬，將數(shù)值模擬的結(jié)果與重建算法得到的結(jié)果進(jìn)行對比分析，重建結(jié)果表明本文的算法能夠較準(zhǔn)確地重建數(shù)值模擬的數(shù)據(jù)，初步證明了本算法的可行性。最后將本算法與基于主成分分析的濃度重建算法的重建誤差進(jìn)行對比，結(jié)果表明了本算法相對于基于主成分分析的濃度重建算法在三維濃度場重建上具有一定的優(yōu)越性。

本文中的重建誤差較大，所使用的測點數(shù)目也較多，后期研究將對提升重建精度和降低測量點數(shù)目進(jìn)行進(jìn)一步探索。

引用

[1] 王雪梅,劉石.基于TDLAS技術(shù)的CO_2濃度測量[J].化工自動化及儀表,2016,43(11):1158-1161.

[2] 丁喜波,陳晨,張任,等.基于超聲波相位差的氣體濃度測量方法[J].高技術(shù)通訊,2014,24(2):189-192.

[3] 李元齊,沈祖炎.本征正交分解法在曲面模型風(fēng)場重構(gòu)中的應(yīng)用[J].同濟(jì)大學(xué)學(xué)報(自然科學(xué)版),2006(1):22-26.

[4] QIN L,LIU S,LONG T,et al.Wind Field Reconstruction Using Dimension-reduction of CFD Data with Experimental Validation[J].Energy,2018,151(May 15):272-288.

[5] SUN S X,LIU S,LIU J,et al.Wind Field Reconstruction Using Inverse Process with Optimal Sensor Placement[J].IEEE Transactions on Sustainable Energy,2018,10(3):1290-1299.

[6] 陳敏鑫,劉石,孫單勛,等.機(jī)器學(xué)習(xí)算法在溫度分布重建中的應(yīng)用[J].自動化儀表,2020,41(1):95-100+105.

[7] 陳敏鑫.基于降維和深度學(xué)習(xí)方法的溫度分布重建[D].北京:華北電力大學(xué)(北京),2021.