黃河下游泥沙長距離管道輸送目標優化方法

劉一鳴 姜晗琳 李振山

摘?要：為緩解目前黃河下游管道輸沙距離短、能效低等問題，針對長距離管道輸沙問題，基于阻力損失與流速的變化特性，建立了以輸沙量偏差和流速偏差為目標函數，以滿足泥沙濃度、流速、管徑、輸沙量等方面要求為約束條件的多目標優化模型。通過黃河下游管道輸沙實例，選用NSGA-Ⅱ遺傳算法進行多目標尋優，得到的Pareto最優解集和最優前沿在取值空間內呈連續均勻分布，且輸沙量偏差向目標函數高效收斂。根據泥沙體積分數不同將所得解集分為3組輸沙方案，以適應動態變化的水沙條件和工程需求。

關鍵詞：泥沙;黃河下游;管道輸送;多目標優化;遺傳算法

中圖分類號：TV67;TV882.1?文獻標志碼：A

doi：10.3969/j.issn.1000-1379.2021.05.006

Abstract： In order to mitigate the problems of short distance and low energy efficiency for the pipeline sediment transport in the lower Yellow River floodplain， a multi-objective optimization model was built for pipeline design for the sediment transport in the Yellow River floodplain. The model consisted of two objective functions （sediment discharge deviation and flow velocity deviation） and four constrains （sediment concentration， flow velocity， pipe diameter and sediment discharge）. NSGA-Ⅱ genetic algorithm was chosen for resolving multi-objective optimization model. The case study shows that the Pareto optimal solution set and the optimal frontier are distributed continuously and uniformly in the value space. Moreover， the deviation function of sediment transport has better convergence. According to the sediment concentration， the solution set was divided into three groups to adapt to the variation of water and sediment conditions and engineering requirements.

Key words： sediment; Lower Yellow River; pipeline transportation; multi-objective optimization; genetic algorithm

黃河下游因泥沙淤積而形成地上懸河，危害甚大[1]。緩解泥沙淤積的有效方法之一是泥沙的綜合利用，泥沙可用于加固堤防、淤灘造地、礦坑填埋、填海造陸等，這些利用方式往往涉及泥沙的輸送。管道輸沙采用多級泵串聯形式，理論上可以實現長距離輸送[2]。近年來管道輸沙研究主要關注某一固定輸沙參數和水沙條件的優化，如計算臨界不淤流速[3]、建立管道輸沙阻力損失[4]等參數模型、輸送參數和動力系統的設計計算等[5]。這些研究在一定程度上解決了區域內管道輸沙設計問題，但仍存在不足。單一目標優化思路雖然可以優化某一輸沙參數，但可能限制了其他輸沙性能參數的優化，而且多基于具體案例的輸沙目標給出設計方案，普適性較差。因此，開展管道輸沙多目標優化方法的系統性研究很有必要。本文將優化思路引入管道輸沙設計中，建立適用于長距離管道輸沙的多目標優化模型，在已有研究的基礎上給出參數的確定方法，并將模型應用于黃河下游輸沙方案優化，根據優化解集的分布規律提出輸沙方案。

1?優化模型

1.1?目標函數

從管道輸沙的實際需求出發，確定優化模型包括兩個目標函數：輸沙量偏差值和流速偏差值。

（1）輸沙量偏差值。為滿足輸沙任務的設計要求，管道應具備足夠的輸沙能力，衡量指標是輸沙量達到目標輸沙量。因此將實際輸沙量與目標輸沙量之差的絕對值最小作為衡量輸沙能力的目標函數，即輸沙量偏差最小，表示為

式中：QΔ為輸沙量偏差，m3/h;Q為實際輸沙量，m3/h;Qt為目標輸沙量，m3/h。

（2）流速偏差值。不同于純流體的輸送，管道輸沙的動力系統做功除了維持流體流動外，還要用于輸送泥沙，包括泥沙與管壁的摩擦損耗以及維持顆粒不沉降和克服顆粒間碰撞造成的機械能損失。鑒于管道內泥漿輸送能量損耗的微觀機理極其復雜，通常將管道阻力損失簡化為管道液體阻力損失和固體附加阻力損失之和[6]，即

式中：im為管道阻力損失;iw為液體阻力損失;is為固體附加阻力損失。

式（2）中的阻力損失im可以通過經驗公式計算，此處選擇Wilson阻力損失模型[7]計算阻力損失：

式中：Cv為泥沙的體積分數;s為泥沙體積質量;u為設計流速;u50為參考流速;Dt為泥沙分布系數;λw為清水阻力系數;D為管道直徑;g為重力加速度，取9.81 m/s2。

研究發現[3，6，8]，管道輸沙的阻力損失與流速的關系：隨流速u增大，阻力損失im與流速的關系曲線呈下凹形，即先減小后增大。當流速達到某個最優流速時，管道阻力損失最小。因此管道輸沙系統應將該流速作為目標流速ut，將管道流速u與目標流速ut之差的絕對值最小作為衡量管道阻力損失最小的目標函數，即

1.2?約束條件

黃河下游實際條件會對輸沙過程產生限制，因此需要用約束條件在優化模型中反映此類限制。約束條件分為管道不淤堵條件和邊界條件兩類，它們共同對流速、輸沙量、管徑等參數進行約束。

（1）管道不淤堵條件約束。泥漿在管道內微觀流動形態復雜且受流速影響，當泥漿流速足夠大時，泥沙顆粒隨水懸浮運動;當泥漿流速減小時，部分泥沙顆粒會淤積在管道底部;當泥漿流速進一步減小時，泥沙顆粒沉淀達到閾值后會出現堵管[8]。考慮到臨界流速與目標流速的關系無法確定，將臨界流速作為管道設計流速的約束條件。

理論上臨界流速作為設計流速的下限可以避免管道出現淤堵，然而，臨界流速往往是通過局部試驗或經驗公式計算得到的，且泥漿中會混有泥團等物質[5]，出于安全考慮，實際設計流速會在臨界流速的基礎上留有一定余量，可以通過增加一個系數n來表達，即

式中：uc為臨界流速。

關于式（5）中臨界流速的研究已有數十年。國外學者如Durand[9]將管道底部開始出現泥沙淤積時的極限淤積流速作為臨界流速，后經Wasp等研究[10-12]進一步得到一系列極限淤積流速的經驗模型。我國學者如費祥俊[13]則將管道內垂向泥漿濃度均勻分布所需的不淤流速作為臨界流速。在常用的臨界流速經驗公式中，Wasp公式在Durand公式的基礎上增加了泥沙粒徑對臨界流速的影響，比較全面反映了各參數對臨界流速的影響，所需參數少，因此此處選擇Wasp模型計算臨界流速，公式為

式中：d50為中值粒徑。

（2）邊界條件約束。綜合考慮管道運行、輸沙能力和水沙條件等實際情況，模型的邊界約束包括管徑約束、泥沙體積分數約束、流速約束和輸沙量約束。水力輸送管的管徑要根據輸沙量、設計流速和泵的性能確定，增大管徑在一定程度上可以減小阻力，但在流量不變的情況下會減小流速，容易造成泥沙淤積，因此管徑D需要給定上下限Dmax、Dmin。泥沙體積分數可通過改變取沙點和取沙方式進行調節，模型中泥沙體積分數Cv的取值范圍要在給定的上下限內，用Cvmax、Cvmin分別表示泥沙體積分數的上、下限。臨界流速作為流速的下限已經在式（6）中考慮。流速越大對泵等動力設備性能要求越高，因此設計流速的上限應在合理范圍內，即u≤umax，其中umax為流速上限。作為衡量輸沙能力的參數，輸沙量遠低于設計目標、無法完成輸沙任務的方案，雖然針對模型具有可行性，但沒有實際意義，因此通過約束輸沙量下限的方式剔除明顯不合理的輸沙方案，約束條件為輸沙量≥Qsmin，其中Qsmin為輸沙量下限。

綜上，將優化方法用于管道輸沙過程，分別將管道輸沙量和流速與目標值的偏差作為目標函數，確定管道運行的邊界約束條件，建立包括模型參數在內的管道輸沙優化模型。

2?參數確定與計算方法

為驗證模型的可行性，用實際算例對模型進行檢驗。先確定模型參數，再選擇多目標優化方法求解模型獲得Pareto最優解的解集，最后對Pareto最優解進行分析并評估輸沙方案。

2.1?參數確定

求解優化模型要先給出基礎參數的取值，管道阻力系數可以通過管道阻力損失計算。計算式（3）中參考流速u50需要伴隨速度W、沉降速度vt、無因次粒徑d*等參數，各參數計算公式見表1，其中：μw為清水黏度;ρw為清水密度;ψ為泥沙形狀系數;ε為絕對粗糙度。

在求解目標流速ut的過程中，由于柯爾布魯克方程需要迭代計算清水的阻力系數λw，因此難以將式（4）求導變換后直接給出ut的表達式，此處取流速步長為0.001 m/s分別試算im，將試算結果中im最小值對應的流速作為該工況下的ut。

以黃河下游調節灘槽地形、清淤疏浚的管道輸沙為例，對上述優化方法進行算例分析。選取黃河下游河床質泥沙作為輸送介質，參考實際水沙、地形條件取值以及常見管道水力輸送系統的設計方法[5，14-17]，定義邊界約束、水沙參數、管道參數等模型參數，結果見表2。需要說明的是：由于工程中的管道均有標準化管徑，因此模型中管徑D簡化為在約束邊界內以0.05 m為步長取離散值;根據泥沙的體積分數、管徑、流速的取值邊界計算出輸沙量的有效取值為3.21～636.17 m3/h，為獲得適用范圍大且輸沙能力較強的設計方案，設置目標輸沙量Qt和輸沙量下限Qsmin分別為500、100 m3/h;為滿足長距離輸沙的要求，取管道長L=3 000 m。

2.2?求解方法

遺傳算法是一種模擬生物進化過程的優化算法，具有良好的適應能力和全局優化能力[18]。非支配排序遺傳算法（NSGA）是遺傳算法面向多目標優化問題的主要應用之一[19]。NSGA-Ⅱ算法改進了NSGA算法過于復雜等缺陷，提高了算法的計算速度和魯棒性，算法新增的精英策略可以擴大模型取值范圍，使得Pareto最優解的分布趨于均勻化，近年來得到廣泛關注，因此將NSGA-Ⅱ算法作為多目標管道輸沙模型的求解算法。

使用modeFRONTIER多目標優化軟件完成求解過程，modeFRONTIER軟件包括一系列多目標優化算法，同時還提供了MATLAB、Excel等計算軟件接口，便于集成優化。優化模型在modeFRONTIER軟件模塊設置中可以分為：決策變量、目標函數、約束條件、模型求解算法4個部分。設置泥沙體積分數Cv、管徑D、設計流速u共3個決策變量，NSGA-Ⅱ算法參數設置：初代種群規模為100，最大遺傳代數為300代，交叉概率為0.9，模型的數值計算由MATLAB完成。

管道輸沙模型優化流程見圖1，可分為多目標優化模型求解、參數分析及方案比選兩步，先用NSGA-Ⅱ算法求解管道輸沙多目標優化模型，若兩目標函數無法同時達到最優解，則將uΔ、QΔ無法再同時改進的解集作為Pareto最優解集，接著對計算結果進行機械能衡算求出動力設備的揚程，最后將Pareto最優解集和揚程匯總分類得到輸沙方案并對方案進行比選。

3?計算結果與討論

3.1?解集空間分布

NSGA-Ⅱ遺傳算法優化后共得到30 000個子代，其中Pareto最優解集由4 446個子代組成。決策變量的Pareto最優解集構成的曲面和目標函數的Pareto最優解集（Pareto最優前沿）空間分布見圖2、圖3。圖2的曲面和圖3的最優前沿呈上凸形，分布較均勻，連續性好，最優解之間不存在互相支配的關系，說明NSGA-Ⅱ算法對于該模型有較強的Pareto解搜索能力。

圖2中所表示的決策變量Cv、D、u的Pareto最優解分別為0.103～0.271、0.45～0.50 m、2.09～2.77 m/s。從Cv的Pareto最優解集取值范圍可以看出，模型的Pareto最優解覆蓋了實際條件下各種不同的泥沙體積分數。從圖3可以看出，兩個目標函數由于優化方向相反，因此模型不存在單一的最優解，輸沙量偏差和流速偏差分別為0.02～374.47 m3/h、0.96～1.22 m/s。分別取解集中的兩個極值，即輸沙量偏差QΔ和流速偏差uΔ的最優個體，計算出對應的設計參數輸沙量Qs和阻力損失im分別為499.98 m3/h、8.3×10-3。受流速邊界、目標流速和臨界流速等條件限制，圖3中uΔ取值范圍較小，可以通過調整流速邊界、現場測定臨界流速等方式來改善目標函數uΔ的取值并提高結果的準確性。

需要指出的是，得到的部分Pareto最優解集和最優前沿（圖2、圖3圈內部分）已經十分逼近目標函數QΔ極值處，算法的收斂性和搜索性能有一定下降，造成了圖2曲面的平滑性和圖3最優前沿的連續性下降。

3.2?解集頻率分布

根據Pareto最優解的頻率分布，研究優化結果向目標函數的收斂情況，圖4、圖5分別為QΔ、uΔ的Pareto最優解頻率分布直方圖，直方圖組距分別取50 m3/h、0.04 m/s。從圖4可以看出，QΔ的Pareto最優解向橫軸左端聚集，符合式（2）規定的QΔ最小化的目標。uΔ的分布較為分散。分別計算uΔ和QΔ的分布特點得出：Pareto最優解在uΔ前25%的頻率為0.19，對應在QΔ前25%的頻率為0.39，可見QΔ向目標函數收斂性更好。相較于流速偏差uΔ，輸沙量偏差QΔ目標函數更具有效性。從子代的分布特性可以看出，在搜索可行解時，種群更易于向達到輸沙量偏差目標的方向演化，說明模型中較多的Pareto最優解集表現為適當提高流速同時增大管道阻力，以此完成輸沙任務。

3.3?輸沙方案比選

根據模型Pareto最優解集和最優前沿分別計算輸沙量Q、流速u、阻力損失im等參數。此外，在實際輸沙過程中，泵揚程He也是值得關注的設計參數，此處使用流體的伯努利方程能量衡算獲得輸沙泵的揚程，管道阻力采用范寧公式計算。

式中：He為揚程;Δz為高差;∑ζ為管道局部阻力系數。

根據黃河下游灘區實際地形取Δz=6 m。管道局部阻力系數∑ζ的計算：管道出入口和管件閥門包括：入口（1個）、出口（1個）、45°彎管（3個）、90°彎管（3個）、截止閥（3個），其對應的局部阻力系數ζi分別為0.50、1.00、0.35、0.75、6.00，對于泥沙漿液，將上述局部阻力系數乘以1.5后加和作為管道局部阻力系數∑ζ[5]。

根據泥沙體積分數不同將模型的Pareto最優解和最優前沿分類得到輸沙方案，并對方案的適用性進行評估。按照泥沙體積分數Cv從小到大的順序，將Pareto最優解集分為3個規劃方案（見表3），從方案整體可以看出，隨著Cv的增大，Q增大，即系統輸沙能力提升，受臨界流速約束的u同時增大，造成管道阻力損失im增大，為了維持管道正常運行，輸沙泵的揚程He也須增大。此外，管徑D較大，說明在該模型的不同輸沙條件下，選擇大管徑有利于輸沙。

對比不同方案可以看出，方案1的管道阻力損失較小，設計流速u較小，但輸沙量也較小，在目標輸沙量Qt取500 m3/h的條件下，完成輸沙量任務低于50%，可見該方案適用于動力條件有限、輸沙規模小的場合。方案3可以完成73%至100%的輸沙目標，但設計流速大，泥沙體積分數大，阻力損失也較大，輸沙泵所需揚程高，必要時需考慮安裝接力泵，適用于輸沙設備動力充足、輸沙目標優先的場合。方案2的管徑D固定為0.50 m，可以完成50%至72%的輸沙任務，設計流速u和泵揚程He取值相對合理。

從表3中不難發現，沒有得到輸沙量Q遠大于Qt的輸沙方案，這可能是目標函數及模型邊界設置造成的，可以通過調整式（1）中QΔ的定義和約束邊界，或者將實際目標輸沙量適當增大作為模型的Qt等方法解決該問題。此外，若實際工況下泥沙體積分數Cv難以調節，則可以先根據當地水沙條件確定Cv，然后在Pareto最優解集中選擇符合Cv取值的設計參數為設計輸沙方案提供參考。

黃河下游管道輸沙多目標優化模型的臨界流速、阻力損失等參數均采用經驗公式，因此具有一定局限性，可以通過水沙試驗等方式獲得更精確的模型參數與邊界條件。

4?結?論

結合黃河下游實際情況，基于管道兩相流的阻力損失、臨界流速等特征參數提出目標函數和約束條件，建立了管道輸沙的優化模型，并通過算例驗證了該模型的可行性。以輸沙量和流速與目標值的偏差為目標函數，以臨界流速和參數邊界為約束條件，建立了符合實際的管道輸沙多目標優化模型;將黃河下游的河床質泥沙作為輸送介質，利用NSGA-Ⅱ遺傳算法求解模型，獲得了連續均勻的Pareto最優解集，說明NSGA-Ⅱ遺傳算法求解優化模型具有可行性，為優化配置輸沙方案提供了基礎;輸沙量偏差Pareto最優解的分布在前25%的頻率為0.39，可見輸沙量具有很強的目標收斂性;根據Cv從小到大的順序，將Pareto最優解集分為3個輸沙方案，方案可用于變化的大體積分數泥沙輸運，決策者可根據工程需求和決策偏好選擇某一目標最優或多目標最優方案。

參考文獻：

[1]?王兆印，林秉南.中國泥沙研究的幾個問題[J].泥沙研究，2003，28（4）：73-81.

[2]?閆大鵬，李德營，周風華，等.煤礦塌陷區利用黃河泥沙進行治理方案研究[J].中國水土保持，2013（9）：34-36.

[3]?鄒履泰，章少強.兩相流管道臨界不淤流速的預估[J].武漢水利電力學院學報，1990，23（增刊1）：139-144.

[4]?李時，秦毅，李國棟，等.小浪底水庫泥沙管道輸送的阻力損失分析[J].水資源與水工程學報，2016，27（5）：150-157.

[5]?衛國峰，王文東，李赴龍，等.泵淤固堤技術參數選擇與動力系統優化[J].人民黃河，2001，23（5）：3-5.

[6]?費祥俊.固體管道水力輸送摩阻損失的預測[J].水利學報，1986，17（12）：20-29.

[7]?WILSON K C，LULE T U， ARKITEKTUR O V， et al. Slurry Transport Using Centrifugal Pumps[M]. US： Springer， 2006： 48-157.

[8]?鄧祥吉.管道輸沙阻力特性研究[D].南京：河海大學，2005：7-8.

[9]?DURAND R. Basic Relationships of the Transportation of Solids in Pipes-Experimental Research[M]. Minnesota： Proc Minnesota Intern Hyd Conv， 1953： 89-103.

[10]?WASP E J， AUDE T C， SEITER R H， et al. Hetero-Homogeneous Solids/Liquids Flow in the Turbulent Regime[C]//Advances in Solid-Liquid Flow in Pipes and its Application， Amsterdam： Elsevier Inc.， 1971：199-210.

[11]?SHOOK C A. Pipelining Solids： the Design of Short-Distance Pipelines[M]. Canada： Soc. Chen. Engin， 1969：13-25.

[12]?瓦斯普.固體物料的漿體管道輸送[M].北京：水利出版社，1980：5-8.

[13]?費祥俊.漿體管道的不淤流速研究[J].煤炭學報，1997，22（5）：86-90.

[14]?閉治躍.挖泥船泥漿管道輸送系統效率優化與控制研究[D].杭州：浙江大學，2008：25-48.

[15]?侯素珍，李婷，郭彥，等.西柳溝河床泥沙組成分析[J].人民黃河，2018，40（7）：7-10，14.

[16]?王志魁，劉麗英，劉偉.化工原理[M].北京：化學工業出版社，2010：12-30.

[17]?和瑞莉，李靜，張石娃.黃河流域泥沙密度試驗研究[J].人民黃河，1999，21（3）：8-10，49.

[18]?王麒，李松青，王心義，等.基于遺傳算法的水源熱泵系統抽灌量優化配置[J].水利學報，2019，50（6）：743-752.

[19]?高媛.非支配排序遺傳算法（NSGA）的研究與應用[D].杭州：浙江大學，2006：1-10.

【責任編輯?呂艷梅】