黃河下游泥沙長距離管道輸送目標優(yōu)化方法

劉一鳴 姜晗琳 李振山

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

關(guān)鍵詞：泥沙;黃河下游;管道輸送;多目標優(yōu)化;遺傳算法

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

1?優(yōu)化模型

1.1?目標函數(shù)

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

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

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

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

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

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

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

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

1.2?約束條件

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

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

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

式中：uc為臨界流速。

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

式中：d50為中值粒徑。

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

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

2?參數(shù)確定與計算方法

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

2.1?參數(shù)確定

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

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

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

2.2?求解方法

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

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

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

3?計算結(jié)果與討論

3.1?解集空間分布

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

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

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

3.2?解集頻率分布

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

3.3?輸沙方案比選

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

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

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

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

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

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

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

4?結(jié)?論

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

參考文獻：

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

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

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

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

[5]?衛(wèi)國峰，王文東，李赴龍，等.泵淤固堤技術(shù)參數(shù)選擇與動力系統(tǒng)優(yōu)化[J].人民黃河，2001，23（5）：3-5.

[6]?費祥俊.固體管道水力輸送摩阻損失的預(yù)測[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]?閉治躍.挖泥船泥漿管道輸送系統(tǒng)效率優(yōu)化與控制研究[D].杭州：浙江大學，2008：25-48.

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

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

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

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

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

【責任編輯?呂艷梅】