基于粒子群算法的彎道環流式流道多目標優化

馬睿佳，魏正英，陳雪麗，馬勝利(西安交通大學機械制造系統工程國家重點實驗室，西安 710049)

迷宮型灌水器具有結構相對簡單、制造方便、成本低，又具有一定的壓力補償性等諸多優點而得到廣泛應[1]。迷宮型灌水器主要由3部分組成：進口柵格、迷宮流道和出口。其中迷宮流道利用復雜的壁面邊界使其中的水流紊亂，從而引起能量損失，水流在流道內的流動近似表現為紊流，是此類灌水器消能的主要部位，它的結構特性在很大程度上決定著灌水器的2大性能：水力性能和抗堵性能。

傳統的求解方法需要對多個目標進行合并,需要選取適當的折中系數,有很強的先驗性,增加了決策者決策的難度，如灰靶理論[2]等。為了更好地解決復雜的多目標優化問題，研究人員引進了多目標進化算法，其中，粒子群算法，也稱粒子群優化算法(Particle Swarm Optimization 簡稱PSO)， 是近年來發展起來的一種新的進化算法。PSO 算法是從隨機解出發，通過多次迭代尋找最優解也稱為非劣解，它也是通過適應度來評價解的品質，它通過追隨當前搜索到的最優值來尋找全局最優[3]。本文通過對兩大性能的定量分析，在正交實驗的基礎上，采用粒子群算法對彎道環流式流道進行多目標優化。

1 彎道環流式流道結構設計

1.1 設計依據

彎道環流現象是指在2個連續且相反的90°圓弧連接處，由于前后兩端壓強的作用剛好相反，會形成2個轉向相反的螺旋狀前進流，在兩個圓弧連接處垂向截面有雙環流現象[4]。流體的能量損失包括局部損失和沿程損失，對于灌水器的流道而言，由于內部流道的長度相對較短，能量損失主要為局部損失。利用彎道環流現象可以增加流動內部的碰撞和漩渦，增加流態內部的紊亂程度，從而有效耗能。借鑒上述現象設計彎道環流式流道如圖1所示。

圖1 彎道環流式流道Fig.1 Bend circulating current emitter 注：1-出口；2-進口；3-流道；4-直線過渡段；5-寬彎段；6-圓弧漸縮段；7-細彎段。

1.2 流道結構參數化

流道單元包括直線過渡段、寬彎段、圓弧漸縮段及細彎段4部分，如圖1所示。寬彎段和細彎段分別為2個同心圓弧構成，且2個部分的圓心位于同一水平上，直線過渡段與2段圓弧外圍進行相同大小的圓弧連接，圓弧漸縮段為過渡輔助。流道單元設計為寬進窄出，在圓弧漸縮段產生漸縮、碰撞及本身的螺旋渦旋損失，進入直線過渡段后又產生漸擴損失。如圖2所示，提取流道中結構參數，寬彎段的內徑為R1，直線過渡段兩端的過渡圓弧半徑都為R2，細彎段的內徑為R3，直線過渡段寬度為H，2段圓弧的外徑保持不變，寬彎段外徑為1.3 mm，細彎段的外徑為1.4 mm，漸縮段為與寬彎段相切的且必須通過A點的圓弧。

圖2 流道單元結構Fig.2 Channel unit of emitter

2 流道結構正交實驗設計

2.1 正交實驗設計

選擇流道參數R1、R2、R3、H共4個因素，不考慮個因素之間的交互作用，加入誤差項，選用5因素4水平L16(45)正交表，試驗因素水平表見表1，正交實驗設計見表2。共16組，采用數值模擬的方法，通過建模-網格劃分-FLUENT計算的方式對流道性能參數進行提取分析，建立流道結構參數與灌水器2大性能之間的目標函數。

表1 正交實驗因素水平 mmTab.1 Factors and levels of orthogonal experiment

表2 L16(45)設計及試驗結果Tab.2 L16(45) design table and experiment results

2.2 數值模型

灌水器流道深度為1 mm,流道進出口采用六面體結構網格劃分，流道部分由于多圓弧部分，采用四面體網格劃分，網格密度為0.1 mm,總網格數為58萬個，細化網格密度分別為0.08和0.06 mm，統計質量流誤差不超過1%，但計算耗時分別增長至2～3倍。流道進口為壓力進口，出口為大氣壓。模型采用標準k-ε模型，近壁面選用標準壁面函數，采用SIMPLE算法，二階離散格式進行離散，收斂精度為10-4。采用離散相模型進行流道內部固液兩相流的研究[4]。

2.3 正交實驗結果分析

取流態指數作為水力性能的衡量指標，粒子通過率作為抗堵性能的指標，進行16組數值模擬，分別采用直觀比較分析、方差分析、線性回歸分析等方法對表2進行分析，并得到結構參數與2個指標的回歸模型。

(1)流態指數分析。通過直觀分析可知各結構對流態參數影響為R1>R3>R2>H，根據流態指數因素-指標圖3，選擇流態指數的極小值，可知各因素的最佳組合為R1=0.6 mm，R2=0.65 mm，R3=0.71 mm，H=0.85 mm，此組合沒有在正交組合中出現。對其進行方差分析，對α=0.1，查表可知F0.1(3,3)=5.39，結果如表3所示，可知R1對結果的影響一般顯著，其他參數對結果影響不顯著。用SPSS軟件對數據進行回歸分析，得到流態指數與各結構參數的回歸模型為：

X=3.022-1.511R1-0.989R2-2.325R3-2.239H+

1.243R21+0.781R22+1.667R23+1.319H2

(1)

表3 流態指數方差分析結果Tab.3 Variance analysis results of flow index

圖3 流態指數因素-指標圖Fig.3 The factor-index chart of flow index

模型的相關系數R=0.965，R2=0.931，線性關系顯著。將通過回歸模型計算的流態指數與正交試驗的結果進行比較，二者的相對誤差在2.6%以內，回歸模型精度較高。

(2)粒子通過率分析。同理，通過直觀分析可知流道結構參數對對粒子通過率的影響順序為H>R1>R2>R3。具有最大的粒子通過率的結構組合為：R1=0.6 mm，R2=0.60 mm，R3=0.74 mm，H=0.85 mm。在相關系數R方為0.878的顯著水平下，回歸分析得到的粒子通過率與流道結構參數的回歸模型為:

K=3.909-3.215R1-0.067R2+3.482R3-8.602H+

3.09R21+0.133R22-2.243R23+5.168H2

(2)

3 基于粒子群算法的流道結構優化

3.1 多目標粒子群算法

粒子群算法是一種隨機的基于種群的優化策略，通過模擬鳥群覓食的過程建立的一種優化算法。粒子群算法中，每只鳥作為一個粒子，每個粒子有位置和速度兩個屬性，在飛行過程中根據自身找到的距離食物最近的位置和整個群體找到最近位置去不斷改變自身的前進方向，最后整個群體都會趨向同一個地方，即食物的位置。在原有的粒子群算法基礎上，Shi[6]等提出了慣性權重粒子群算法，也稱標準的粒子群算法，包括了3個部分，粒子自身的慣性部分、自身認知部分和社會部分。該算法的速度更新公式如下：

vid(t+1)=ωvid(t)+c1ri1,t[Xpbid(t)-xid(t)]+

c2ri2,t[Xgbid(t)-xid(t)]

(3)

式中：ω為慣性權重，描述粒子慣性對當前速度的影響；c1、c2為加速系數；ri1,t、ri2,t為0到1之間的隨機數。位置更新公式如下：

xid(t+1)=xid(t)+vid(t+1)

(4)

之后，Clerc[7]等提出了一種帶有收縮因子的粒子群算法，在保證算法收斂的同時提高算法的收斂速度。本文采用這種算法對速度進行更新，該算法的描述如下：

vid(t+1)=χ{vid(t)+c1ri1,t[Xpbid(t)-xid(t)]+

c2ri2,t[Xgbid(t)-xid(t)]}

計算取c1=c2=2.05，則φ=4.1，χ=0.729。由式(3)和式(5)可知，2種算法從表達形式來說基本一致，式(5)是式(3)的變形。

多目標粒子群算法(multiobjective particle swarm optimization,MOPSO)中以2004年Coello[8]提出基于精英級策略的MOPSO算法最具有里程碑意義，本文采用該算法進行灌水器結構的優化計算。該算法的具體步驟如下。

(1)隨機初始化種群，計算各個粒子的適應度值，將其中的非劣解及其位置信息存入外部檔案庫REP。

(2)確定粒子初始的個體最優Xpbi(0)和全局最優Xgbi(0)。

(3)迭代次數增加，每迭代一次，更新粒子的速度和位置，其中個體最優位置通過比較單個個體經歷過的歷史最優位置，同時為每一個粒子選取全局最優，全局最優采用自適應網格法和輪盤選擇法[9]從外部檔案庫REP中選取一個。

(4)迭代次數增加，直至滿足終止條件。最終外部檔案庫REP保存的就是所需要的最優解。

外部檔案庫REP的作用機制包括以下幾個方面：首先，每一次迭代得到的全部粒子中的非劣解與上一次迭代得到的REP內的值進行比較，若新的非劣解被REP內原有的任一解支配，則無法進入REP內；若REP內的所有值都無法支配新的非劣解，則新的非劣解進入；若新的非劣解可以支配其中的某個非劣解，則原來的非劣解丟棄，新的非劣解進入REP內。當外部檔案庫REP達到數量上限時，采取二次準則對里面的元素進行保留，即網格內稀疏位置的元素具有優先保留權。

3.2 基于MOPSO的彎道環流式流道結構優化

灌水器流道多目標優化的目標函數由上文得到如下：

MinX=3.022-1.511R1-0.989R2-2.325R3-

2.239H+1.243R21+0.781R22+1.667R23+1.319H2

MaxK=3.909-3.215R1-0.067R2+3.482R3-

8.602H+3.09R21+0.133R22-2.243R23+5.168H2

其中，結構參數的約束范圍為：

(6)

由于優化目標是找到具有優良的水力性能及抗堵塞性能的小流量灌水器，故將流量作為額外的約束標準，流量回歸方程由上述正交實驗獲得：

Q=2.319-0.803R1-0.253R2-0.343R3+0.354H

1.5

(7)

其中設置種群規模為100個，迭代次數為100次，分別比較經過初始代、30次、70次、100次迭代后的非劣解集如圖4所示，可知粒子通過率K隨著流態指數x的增大而增大，基本上保持為線性關系。而理想的結果為粒子通過率越大且流態指數越小越好，故可知兩者是相互抵制的，在得到非劣解集后，仍需要人為選擇兩者關系的權重來確定相應的組合。分別在迭代第100次得到非劣解集中選取5個點，如圖5所示，并列出對應的結構參數的組合，如表4所示。

由表4可知，并對其中非劣解集中的點的結構參數及目標函數值分析可知，R1的大小基本為0.6，其他3個參數的可浮動量的大小分別為H最大，R2次之，R3最小。說明R1的取值對2大性能的影響很大，只有在取值為0.6的時候才可以得到最優值，而非劣解集的前沿分布主要由H和R2的大小決定。同時，對比前面的正交實驗結果可知，前面的最佳組合基本圍繞在參數變量的附近，可知正交分析與多目標粒子群算法的正確性。得到非劣解集后，應當根據設計需求 期、抽穗期、灌漿期的適宜灌溉控制下限標準(占田持%)分別為45%、65%、70%、65%，且在生產中盡量避免出現長期連續干旱現象。

圖4 迭代過程中非劣解集的進化過程Fig.4 Evolutionary process of non dominated solution set in the iterative process

圖5 取點位置圖Fig.5 Location of point

表4 非劣解參數信息Tab.4 Parameters of non dominated solution

□

[1] 路振廣，孟春紅，王曉麗.基于遺傳算法的冬小麥優化灌溉制度研究[J].人民黃河，2009,31(11)：76-78.

[2] 劉 鈺，汪 林，倪廣恒，等.中國主要作物灌溉需水量空間分布特征[J].農業工程學報，2009,25(12)：6-12.

[3] 高 陽，黃 玲，李新強，等.開花后水分脅迫對冬小麥旗葉光合作用和保護酶活性的影響[J].水土保持學報，2013，27(4)：201-206.

[4] 申孝軍，孫景生，劉祖貴，等.灌水控制下限對冬小麥產量和品質的影響[J].農業工程學報，2011，26(12)：58-65.

[5] 徐建文，梅旭榮，居 輝，等.黃淮海地區冬小麥關鍵生育期不同灌溉水平對產量影響的模擬[J].作物學報，2014，40(8)：1 485-1 492.

[6] 姚 寧，宋利兵，劉 健，等.不同生長階段水分脅迫對旱區冬小麥生長發育和產量的影響[J].中國農業科學，2015，48(12)：2 379-2 389.

[7] 王俊儒，李生秀.不同生育時期水分有限虧缺對冬小麥產量及其構成因素的影響[J].西北植物學報，2000，20(2)：193-200.

[8] 閆永鑾，郝衛平，梅旭榮，等.拔節期水分脅迫-復水對冬小麥干物質積累和水分利用效率的影響[J].中國農業氣象，2011，32(2)：190-195.

[9] 霍治國，李世奎，白月明，等.冬小麥中輕度水分脅迫的增產節水效應研究[J].自然資源學報，2003，18(1)：58-66.

[10] 銀敏華，李援農，周昌明，等.調虧灌水和分蘗干擾對冬小麥生長的補償效應[J].應用生態學報，2015，26(10)：3 011-3 019.