基于MATLAB與COMSOL聯(lián)合仿真的梯形迷宮滴頭流道優(yōu)化

胡宇祥，彭軍志，殷 飛，劉喜峰，李 娜

（吉林農(nóng)業(yè)科技學院，吉林 132101）

0 引 言

滴灌灌水器是滴灌系統(tǒng)中最為關鍵的元件，其使用過程中蒸發(fā)損失小，不產(chǎn)生地面徑流，幾乎沒有深層滲漏。隨著農(nóng)業(yè)用水形勢逐年嚴峻以及國家“節(jié)水增糧行動”項目的大力推廣，滴灌技術成為最有效的灌溉方式，是保障增產(chǎn)增收的重要手段[1]。滴灌灌水器結構尺寸對滴灌系統(tǒng)灌水的均勻性、抗堵塞能力影響顯著。過大的流道結構參數(shù)將造成灌溉均勻度下降，過小的流道結構參數(shù)將造成灌水器流道堵塞，嚴重時會使整個系統(tǒng)無法正常工作[2-3]。因此，研究流道結構參數(shù)與滴頭水力性能之間的關系一直是節(jié)水灌溉領域的研究熱點。

目前，國內(nèi)外學者已在滴灌灌水器流道結構參數(shù)優(yōu)化方面取得了一系列的研究成果，顏廷熠等利用正交試驗法，研究了流道結構參數(shù)對滴頭水力性能的影響[4-6]；王文娥利用 FLUENT 軟件對內(nèi)鑲齒形片狀滴頭的水力性能進行了數(shù)值模擬，獲得了滴頭壓力-流量關系式并分析了滴頭內(nèi)速度、壓力分布[7]；李永欣等建立了迷宮滴頭的CFD數(shù)值模型，并對滴頭的壓力流量關系、流道內(nèi)部的壓力和流速分布進行了數(shù)值模擬計算，利用原型滴頭和滴頭放大模型實測值對模型和模擬計算結果進行了試驗驗證[8]；陳愛平利用 GAMBIT軟件建立了輸送管的仿真模型，通過 FLUENT 軟件的數(shù)值模擬計算，得到了不同滴灌裝置輸送管結構參數(shù)下的流場分布特性[9]；Souza等通過試驗方式，設計了一種灌溉滴頭，測定了滴頭流量系數(shù)和流態(tài)指數(shù)[10]；Li等利用多元回歸分析，建立了流量與流道幾何參數(shù)之間的線性回歸方程[11]；Wei等利用CFD對滴頭水利特性進行了仿真研究[12]。目前，上述研究多是利用正交試驗方法，將流道結構參數(shù)設計成多因素多水平的試驗，通過CFD 方法或者試驗法獲得滴頭的宏觀水力特性，即壓力與流量間的關系，分析流道結構參數(shù)與流態(tài)指數(shù)之間關系[13]。但是通過正交試驗設計進行數(shù)據(jù)分析后得到的優(yōu)選結果，只是因素間不同水平的相互組合，難于確定數(shù)據(jù)變化規(guī)律。水平設計的少，則自變量定義域不完整，試驗設計具有隨機性，結果不夠精確[14-15]；如果水平設計的多，則試驗數(shù)據(jù)體量大，試驗次數(shù)多，尤其是費用昂貴的試驗，弊端更加凸顯。因此，需要一種能夠快速準確求解灌水器流道尺寸的計算方法。

鑒于此，本文提出聯(lián)合MATLAB與COMSOL仿真軟件進行滴灌灌水器流道參數(shù)優(yōu)化設計，采用COMSOL軟件對梯形迷宮流道灌水器進行仿真，求解最佳流態(tài)指數(shù)下的滴灌灌水器結構參數(shù)，研究流道齒參差值、齒高、單元數(shù)、流道轉角、齒間距這5個變量對灌水器流道參數(shù)的影響，并在此基礎上建立數(shù)學預測模型，以期為快速確定灌水器流道結構參數(shù)提供理論依據(jù)以及科學參考。

1 研究方法

1.1 滴灌灌水器流道數(shù)學模型

1.1.1 設計變量

迷宮形滴頭相比貼片形式滴頭抗堵性能好，灌水均勻，滴頭整體性強。Netafim公司作為滴灌技術的發(fā)明者，一直在引領滴灌領域技術變革。本文以該公司迷宮滴頭流道參數(shù)的取值范圍作為參考依據(jù)，其流道結構如圖1所示。

圖1 灌水器流道結構及參數(shù)示意圖Fig.1 Schematic and parameter diagram of emitter channel

本文選取梯形迷宮流道滴頭進行研究，影響其水力性能的關鍵結構參數(shù)包括流道寬度w、流道轉角θ、齒高h、齒間距l(xiāng)、齒參差值j，灌水器流道單元數(shù)n。θ、h、l、j、n這5個參數(shù)可以確定梯形齒道單元平面形狀，流道寬度w按式（1）計算[16-17]。

算例中，流道轉角θ、齒高h、齒間距l(xiāng)、齒參差值j和流道單元數(shù)n是設計變量。根據(jù)已有的研究文獻，各變量取值范圍如表1 所示[13]。

表1 灌水器的結構參數(shù)取值范圍Table 1 Range of structural parameters of emitter

1.1.2 目標函數(shù)

滴灌滴頭的流態(tài)指數(shù)x表示流量對壓力發(fā)生改變的敏感程度，流態(tài)指數(shù)越小，滴頭的水力性能越好[18]。x=1表示流量變化與壓力以相同的百分率變化，x=0表示流量不隨壓力變化而變化。壓力補償式灌水器，流態(tài)指數(shù)趨近于0；非壓力補償式灌水器，流態(tài)指數(shù)介于0～1之間。非壓力補償式灌水器x值越趨近于0.5，流量對壓力敏感程度就越小，水力性能越優(yōu)異[19]。因此，建立目標函數(shù)如式（2）所示，要求非壓力補償式灌水器流態(tài)指數(shù)趨近于0.5，保證灌水器水力性能，提高出流均勻性。

式中F為目標函數(shù)；x為流態(tài)指數(shù)，按照SL/T67.1—94要求，可由流量與壓力關系[20]計算得出：

式中Q為滴頭流量，L/h；H為工作壓力， kPa；x為流態(tài)指數(shù)；Kd為流量系數(shù)。從0.5～1.5倍工作壓力（100 kPa）范圍內(nèi)，均勻選取 9個工作壓力值作為自變量，輸入COMSOL軟件作為滴頭進口壓力，通過仿真計算出滴頭出口流量Q作為因變量，利用MATLAB進行最小二乘法計算，建立基于乘冪函數(shù)形式的流量與壓力關系，得到回歸系數(shù)Kd和x。

1.1.3 邊界條件

滴頭流道進水口設為工作壓力，出口設為當?shù)卮髿鈮海嬎愠隹诹魉佟３跏紶顟B(tài)設置為滴頭內(nèi)充滿水，滴頭以外為空氣，空氣壓力為大氣壓[7]。本文采用精度高、計量小的標準壁面函數(shù)，減小梯形滴頭流道的壁面對紊流的影響[21-22]。

1.2 遺傳算法原理

與傳統(tǒng)的優(yōu)化算法相比，遺傳算法生成初始群體開始搜索，按照基因的表達方式，目標函數(shù)值轉換得到適應數(shù)值的信息，而不需要借助其他輔助信息以及參數(shù)。遺傳算法有 3個基本操作：選擇算子、交叉算子和變異算子，各個算子之間即相互配合又相互競爭[23]。遺傳算法作為一種全局優(yōu)化算法，相對于粒子群算法、蟻群算法和魚群算法等群智能算法，其基本原理簡單易懂，初始條件少，廣泛應用于不同領域。本文利用MATLAB自帶的遺傳算法工具箱進行求解。

1.3 COMSOL與MATLAB聯(lián)合仿真在滴灌灌水器流道優(yōu)化中應用

COMSOL（COMSOL Multiphysics）是多物理場仿真軟件，可用于溫度、電磁、流體等多物理場之間的相互耦合用。 COMSOL中包含livelink for MATLAB模塊，能夠實現(xiàn)與MATLAB軟件之間數(shù)據(jù)的相互傳遞、相互交流，這樣可以利用 MATLAB 及其工具箱進行預處理、模型設定和后處理，進行多物理場模型仿真和遺傳算法優(yōu)化[24]。

COMSOL與MATLAB聯(lián)合仿真優(yōu)化滴頭流道結構流程是，首先運行MATLAB遺傳算法，生成設計變量、設定遺傳優(yōu)化的主要參數(shù)，然后在 COMSOL中完成滴頭流道結構的參數(shù)化建模與計算，最后使用 MATLAB讀取 COMSOL仿真結果，計算目標函數(shù)，完成遺傳算法尋優(yōu)化。COMSOL與MATLAB兩者之間流轉交換步驟如下：

1）運行 MATLAB遺傳算法，生成流道轉角θ、齒高h、齒間距l(xiāng)、齒參差值j和流道單元數(shù)n設計變量。

2）將設計變量和已知變量輸入COMSOL，進行滴頭流道結構的參數(shù)化建模與仿真計算。

3）將仿真計算結果傳給MATLAB，擬合出Q～H公式。

4）設定遺傳算法中選擇、交叉、變異參數(shù)，開始計算目標函數(shù)值。

5）判斷是否滿足收斂條件，如果是，轉向步驟6）；否則，轉向步驟1），生成下一代設計變量。

6）求出最優(yōu)解，計算結束。

1.4 線性回歸分析原理

回歸分析（Regression analysis）是確定 2種或 2種以上變量間相互依賴定量關系的一種統(tǒng)計分析方法。線性回歸分析是最為常用的方法之一，其目的是在因變量和一個或多個自變量之間建立一種關系。該方法可用于因素分析，統(tǒng)計預測，調(diào)整因素等，通過求解標準化回歸系數(shù)，可以分析自變量對因變量的影響程度[25]。求解過程包括多重共線性檢驗和計算標準化回歸系數(shù)。

方差膨脹系數(shù)（Variance Inflation Factor，VIF）是衡量多元線性回歸模型中多重共線性嚴重程度的一種度量。它表示回歸系數(shù)估計量的方差與假設自變量間不線性相關時方差相比的比值。當出現(xiàn)嚴重共線性問題時，會導致分析結果不準確，因而需要及時進行處理。VIF計算如式（4）所示。

式中VIF為方差膨脹系數(shù)；Ri為自變量Xi對其余自變量作回歸分析的負相關系數(shù)。VIF＞5說明存在共線性問題，VIF＞10說明存在嚴重的多重共線性問題，模型構建效果不理想，需要進行處理。

當分析通過多重共線性檢驗后，需要計算標準化回歸系數(shù)，用于比較自變量對因變量的影響程度。標準化回歸系數(shù)絕對值越大說明該自變量對因變量的影響越大。若回歸方程的形式如式（5）所示。

式中Y是估計值；m為變量個數(shù)，取5；多項式系數(shù)bj是非標準化系數(shù)，通過最小二乘法求得。標準化回歸系數(shù)為

式中Beta為標準化回歸系數(shù)；σj為Xj的標準差；σY為Y的標準差。

1.5 數(shù)值模擬

本文利用 COMSOL軟件建立梯形迷宮滴頭流道模型，基于標準k-ε湍流模型進行仿真，計算中對梯形滴頭流道劃分網(wǎng)格。梯形流道結構壁薄且有較多轉角，在對流道劃分網(wǎng)格時，這些結構對網(wǎng)格質(zhì)量的影響十分顯著。為了避免轉角結構給網(wǎng)格劃分帶來困難，采用平面計算模型來研究流道的水力特性，進行平面模型網(wǎng)格劃分[22]。采用的網(wǎng)格大小為0.01～l/20 mm（l為齒間距），因為過密劃分網(wǎng)格對于提高計算精度并不明顯，還會增加計算時間。本文采用漸變方式劃分網(wǎng)格，對于規(guī)則、流動變化緩慢的區(qū)域，采用較稀疏的網(wǎng)格；對于不規(guī)則、流動變化大的區(qū)域采用細密網(wǎng)格，疏密網(wǎng)格間平滑過渡，滴頭進口設置為100 kPa工作壓力。

圖2 迷宮流道的網(wǎng)格劃分Fig.2 Grid generation of labyrinth channel

1.6 試驗方法

為了驗證聯(lián)合仿真計算方法的可行性，評價梯形迷宮滴頭流道灌水器實際的工作性能，利用CNC精雕機制作亞克力材料滴頭測試樣本。按照SL/T67.1—94規(guī)范要求，在0.5～1.5倍工作壓力（100 kPa）范圍內(nèi)，排列25個滴頭按流量由小到大排列編號，取第 3、12、13、23號 4個灌水器為試樣[26-27]。由小到大調(diào)節(jié)滴頭的工作壓力，在測試壓力范圍內(nèi)均勻分布 9個壓力點，用量筒測量4個灌水器每一個壓力點的出水量，時間不少于3 min，試驗重復 3 次，3次所測水量之差不得大于2%，取平均值計算流量，進而得到試樣灌水器的壓力與流量關系曲線。將試樣灌水器每一組壓力下的平均流量進行回歸，繪制灌水器壓力-流量關系曲線。

2 結果與討論

2.1 優(yōu)化結果與驗證

在COMSOL建模基礎上，利用MATLAB遺傳算法工具箱進行優(yōu)化，開展梯形迷宮滴頭流道優(yōu)化聯(lián)合仿真計算。遺傳算法參數(shù)選擇如下：初始種群規(guī)模為 50，遺傳終止代數(shù)MAXGEN= 50，選擇算子采用賭輪盤方式，選擇概率為1/3；交叉算子采用兩點交叉算子，交叉概率為 0.75；變異算子采用高斯變異，變異概率為 0.02。灌水器的結構參數(shù)優(yōu)化結果如表2所示。

表2 灌水器結構參數(shù)優(yōu)化結果Table 2 Optimization results of emitter structure parameters

根據(jù)表2，通過回歸分析得到灌水器流量與壓力關系為Q=0.466 9H0.4999，相關系數(shù)R2為0.99，擬合較好。根據(jù)試驗結果繪制壓力-流量關系曲線如圖3所示，由實測數(shù)據(jù)得到滴頭試樣流量與壓力的函數(shù)關系為Q=0.332 3H0.5576，R2= 0.968。對比仿真優(yōu)化結果，在常壓工作壓力下，實測平均流量普遍偏小。利用SPSS軟件，對仿真優(yōu)化結果與試樣實測結果進行相關性分析，相關性系數(shù)為 0.994，根據(jù)顯著性檢驗方法計算得到P值為0.000 1，小于 0.001，表明結果有較強相關性。根據(jù)GB/T 17187—2009規(guī)定，滴頭流量模擬計算值與實測值之間的平均誤差為6.1%，小于7%，設計的滴頭滿足國標要求[28-29]。

圖3 優(yōu)化結果與實測結果壓力-流量曲線對比Fig.3 Comparison of pressure - flow curve between optimized results and measured results

2.2 滴頭流道流速

在100 kPa工作條件下，優(yōu)化參數(shù)的梯形迷宮流道灌水器流道的流速矢量分布和速度局部放大圖如圖4所示。從滴頭流道流速分布圖可知，各單元之間流速分布一致。根據(jù)滴頭流道速度局部放大圖，可以把梯形滴頭流道內(nèi)部分為3個區(qū)域：低流速區(qū)域A，處于齒道轉角處，速度0.1～1.2 m/s；高流速區(qū)域B，處于軸心位置，速度2.8～3.8 m/s；中流速區(qū)域C，處于A、B區(qū)域之間的過渡區(qū)域，速度1.5～2.5 m/s。由低流速區(qū)域形狀可知，梯形齒尖已經(jīng)趨近于三角形。基于梯形迷宮滴頭流道建立優(yōu)化模型，得到趨于三角形迷宮滴頭流道優(yōu)化結果，可見三角形滴頭相比梯形滴頭優(yōu)勢明顯[18,30]。

圖4 滴頭流道流速分布Fig.4 Velocity distribution in emitter channel

從流道內(nèi)的流速分布圖可知，流道進口段流速較低，是長流道滴頭容易堵塞位置，與試驗研究較為符合[31]。由圖4可知，齒道轉角部位A區(qū)域流速近似為0，該區(qū)域流速慢，顆粒易發(fā)生沉積，從定性角度來說，更易發(fā)生堵塞。以單個齒道為例，A區(qū)域面積占截面面積75.1%，從灌水器流道的結構設計出發(fā)，下一步可將滴頭流道水力性能與抗堵性能作為優(yōu)化目標，將流道形式設計改為圓弧形邊壁，減小低流速區(qū)域面積，從而提高灌水器的抗堵塞性能，優(yōu)化流道設計。

2.3 滴頭流道壓力

最優(yōu)方案的梯形迷宮流道灌水器流道壓力分布如圖5所示，由于整個流道為全封閉狀態(tài)，沒有與外界的能量交換，流道內(nèi)部流體流動過程中產(chǎn)生了沿程水頭損失和局部水頭損失，壓力由進口到出口逐級減小，從100 kPa下降為 0，13個單元中每個單元段水頭損失大約為7.7 kPa，壓降基本相同，壓力沿流道長度呈線性變化。

圖5 滴頭流道壓力分布Fig.5 Pressure distribution in emitter channel

2.4 流道結構參數(shù)對滴頭水力性能的影響

為分析流道轉角、齒高、齒間距、單元數(shù)、齒參差值5個設計變量對流態(tài)指數(shù)x的影響，以仿真優(yōu)化過程中歷代遺傳算法計算結果數(shù)據(jù)為基礎，利用SPSS軟件進行線性回歸分析，結果如表3所示。

表3 流道結構參數(shù)逐步回歸分析結果Table 3 Results of stepwise regression analysis of flow channel structure parameters

由表3得到5個變量與流態(tài)指數(shù)之間的多元線性回歸數(shù)學模型為

從表3可知，最大的VIF為2.531＜5，說明5個變量不存在共線性問題，說明此模型是可靠的。該線形模型的相關系數(shù)R2=0.914，顯著性P值為0.001，說明5個變量與流態(tài)指數(shù)之間線性關系很顯著。通過比較5個設計變量標準化回歸系數(shù)，得到設計變量對流態(tài)指數(shù)x影響程度由大到小的排序為：齒參差值、齒高、單元數(shù)、流道轉角、齒間距。齒參差值、齒高 2個參數(shù)對流態(tài)系數(shù)大小起主要影響作用。齒參差值、流道轉角、齒間距 3個參數(shù)標準化回歸系數(shù)（Beta）為正值，說明增加齒道的齒參差值、流道轉角大小可以增大流態(tài)指數(shù)；齒高、單元數(shù)2個參數(shù)標準化回歸系數(shù)（Beta）為負值，說明減少齒高、單元數(shù)可以增大流態(tài)指數(shù)。由分析結果可知，100 kPa工作條件下，流道結構參數(shù)中的齒參差值對滴頭的水力性能影響較大。

3 結 論

本文通過MATLAB與COMSOL聯(lián)合仿真計算求解出常壓下最優(yōu)梯形迷宮流道灌水器流道結構參數(shù)，流道內(nèi)低流速區(qū)域流速達到0.1～1.2 m/s，過流面積占流道截面的 75.1%，有堵塞風險；高流速區(qū)域流速達到 2.8～3.8 m/s，中流速區(qū)域流速達到1.5～2.5 m/s，流道內(nèi)壓力沿流道長度呈線性變化；仿真結果與實測結果的流量誤差為6.1%；得到5個自變量對因變量影響程度由大到小的排序為：齒參差值、齒高、單元數(shù)、流道轉角、齒間距，建立了流態(tài)指數(shù)與 5個參數(shù)之間的多元線性回歸數(shù)學模型，具有較好的代表性。

本文研究了常壓（100 kPa）工作條件下的流道結構尺寸，未涉及低壓、高壓條件，研究結論有一定的局限性。以流態(tài)指數(shù)作為判斷滴頭水力性能的指標外，還需要對滴頭抗堵塞性能做進一步研究。受滴灌系統(tǒng)工程造價、作物品種等因素影響，應該開展不同工作壓力、不同流量條件下，灌水器流道結構參數(shù)優(yōu)化研究，為開發(fā)新產(chǎn)品、推廣應用提供科學依據(jù)。