自激振蕩脈沖射流減阻特性的多目標優化設計

汪朝暉，馮亞楠，高全杰

（武漢科技大學機械自動化學院，湖北 武漢 430081）

1 引言

高速運動的流體會在壁面邊界層產生強烈的摩擦阻力，是高速流體能量損耗的主要原因。研究表明：高速流射流形態對流道壁面剪切力有較大影響，并且調整自激振蕩腔室結構參數能夠改變射流在邊界層中的擬序結構，可實現流體減阻效果［1-2］。

利用壁面振動引起的湍流減阻效應，控制射流形態，可有效減少流體在壁面處的阻力。文獻［3］利用雷諾數在（9700～19200）的射流，同時對固體壁面進行激振研究，實現了（10～15）%的減阻效果。然后利用儀器對該壁面展翼方向振動時的湍流邊界層變化進行詳細探討。結果表明：射流流體近壁區部分受到壁面振動的影響，其平均速度梯度降低，壁面平均剪切應力下降，從而發現壁面的振動可以有效改善壁面阻力。文獻［4］利用DNS 數值研究發現，流向渦旋具備循環再生的現象。旋渦結構從主射流中獲得能量，這一能量通過旋渦形成了近壁條紋波，同時受到近壁條紋波的影響，又會誘發產生新的旋渦結構。因此，他們認為形成湍流壁面摩擦阻力的關鍵因素是射流中的流向渦旋，而抑制射流中流向渦旋的產生就可以到達降阻的效果。文獻［5］在研究柔性涂層減阻時，將壁面模型設置為多個彈性體，并建立了壁面隨流體波振的環形系統，從而發現振動頻率、射流流速對減阻效果的影響，其減阻率隨著振動頻率和射流流速的提高而變大。文獻［6］人研究了簡諧振動的壁面下的流體運動，利用直接數值模擬，將阻力變化周期分為三個區域進行討論，提出振動引起的負展向渦的產生以及條帶滑移加寬和流向渦減弱是壁面振動產生減阻的主要原因。高速流體管道內沿程阻力主要來源于近壁湍流邊界層，邊界層中的湍流擬序結構會導致管道內的流阻大幅增加，而這種擬序結構主要包括：條帶結構和流向渦結構，條帶結構的向上抬升和向下運動均會導致流阻的急劇增加，而流向渦結構的周期性破碎是導致條帶結構抬升、下降的關鍵原因，因此，抑制湍流擬序結構的周期猝發現象是提高管道內射流減阻效果的關鍵因素。

自激振蕩腔室在一定條件下能夠產生脈沖射流，形成周期性的湍流擬序結構［7］。有效調整射流形態，抑制湍流結構的周期猝發現象，可以實現減阻。文獻［8］研究了自激振蕩腔室的壓力和速度的時均分布特征，并依據這些特征劃分了腔室流域，提出了自激振蕩腔室的脈沖射流形態與自激振蕩腔室結構參數密切相關。文獻［9］通過建立LES模型對自激振蕩腔室的流動特性進行研究，并設計實驗驗證了LES模型的準確性，對比了不同振幅下的射流脈沖性能，發現了自激振蕩腔的腔室長徑比對自激振蕩射流形態有著重要影響。因此，利用自激振蕩脈沖裝置特殊的流體及壁面自振特性，將自激振蕩脈沖理論和壁面振動減阻理論結合起來，通過優化自激振蕩腔室的關鍵結構參數，來抑制湍流流體近壁區的擬序結構的周期猝發現象，能夠有效減少管道流阻。

為了獲得減阻效果最佳的自激振蕩腔室流道結構參數，需要設計有效的多目標優化方法。文獻［10］采用正交試驗法建立初始樣本點，由多目標正交試驗矩陣分析得出各因素權重，通過BP神經網絡與遺傳算法進行協同優化，得出了最優結構參數。文獻［11］采用拉丁式立方法設計實驗生成初始樣本點，近似響應面法擬合出近似模型，運用NSGA-II多目標優化，求解出Pareto最優解集，實現了對最小壓降和最大換熱的多目標優化；文獻［12］采用田口設計，Kriging響應面法構建模型，NSGA-II對板翅式換熱器組進行多目標優化，討論了四因素的二階交互作用對換熱的影響。優化效果關鍵在于兩點：近似模型和優化算法。由于自激振蕩腔室的設計變量與射流減阻特性之間沒有直接精確的數學關系式，CFD模擬與響應面擬合精度取決于計算成本，同時優化算法存在一定計算誤差，加之優化結構的誤差隨優化過程不斷累計。因此，探索更高優化精度的優化方法，研究針對實際問題的有效多目標優化方法極其關鍵。

鑒于此，為了行之有效的優化自激振蕩腔室結構，達到自激振蕩脈沖射流的最佳減阻效果。這里以自激振蕩腔室脈沖射流為研究對象，采用大渦模擬的方法分析了不同結構參數下的自激振蕩腔室在氣體管道運輸中的減阻效果，并通過面心復合設計（FCCD）方法與響應面法結合來建立近似模型，采用NSGA-II算法對自激振蕩腔室無量綱結構參數進行多目標優化，得到最優自激振蕩腔室結構，為自激振蕩脈沖減阻增輸裝置設計提供理論基礎。

2 自激振蕩腔室助推渦減阻效應

流體在自激振蕩腔室出口內高速運動的過程中，流體阻力的形成主要可歸納為兩部分，一部分是流體與固體壁面之間的摩擦阻力，另一部分是流體之間具有相對運動產生的剪切應力。這些阻力會導致流體的自身動能產生耗散，而耗散的動能轉化為流體內的熱能。流體的粘性和慣性和固體不光滑表面的阻礙作用是導致這一現象發生的重要原因。因此，在分析流體流動阻力時，必須對流體的剪切應力以及壁面特性進行著重分析。流體在沿壁面流動時，靠近壁面的粘性底層處流體的速度為零。流體流速沿著垂直方向從零開始逐步增加，然后逐步減緩直至基本不變。故流體在管路中流動可以分為邊界層區域和主流區域。其中，流動主要阻力都集中在邊界層區域，該區域速度變化較大，而主流區域離壁面較遠，速度基本上不變的區域，其流動阻力很小可以忽略。自激振蕩腔室出流管道中單位面積上的剪切應力為［13］：

式中：τ—單位面積上的剪切應力；μ—內摩擦系數（黏滯系數）；du∕dy—速度梯度。其中，速度梯度表示速度沿垂直于速度方向y的變化率，同時它也是直角變形速度，稱為剪切變形速度。因此，式（1）可以理解為：剪切應力與速度梯度成正比。剪切速度應變率可以用來表征剪切應力，同時可以根據模擬計算得到壁面剪切應力推導計算壁面的摩擦阻力［13］：

式中：τ—單位面積上的剪切應力；d—管道直徑；l—管道長度。

在管道直徑和管道長度一定的條件下，壁面摩擦阻力和壁面剪切應力呈正比關系，因此，該摩擦阻力也和速度梯度呈正比關系，這里利用速度梯度變化來表征流體之間的內摩擦力變化。

在經過自激振蕩腔室后，這種基本穩定的速度分布發生改變。自激振蕩腔室出流管道水平射流剪切層渦環，如圖1所示。靠近壁面的反向助推渦是由剪切層分離后形成的。這些反向助推渦夾帶下的高速主射流同固體壁面沒有直接相接觸，起到了類似空氣軸承的作用，將高速流體同固體壁面進行分層。

圖1 自激振蕩室出流管道水平射流剪切層渦環Fig.1 Self-Excited Oscillation Chamber Outflow Pipe Horizontal Jet Shear Layer Vortex Ring

3 數值計算

3.1 控制方程

與其他湍流模型相比，大渦模擬能夠更加精準的捕捉瞬態流動，因此可以采用大渦模擬有效求解自激振蕩腔室周期性瞬時的剪切應力和速度梯度變化。大渦模擬采用濾波函數將湍流瞬時運動中尺度比濾波函數小的渦環濾掉，從而得到大渦流動的運動方程。在非定N-S方程的基礎上進行Favre濾波后得到可壓縮理想氣體的連續性方程、動量方程、能量方程［7］。

3.2 有限元模型

3.2.1 物理模型及關鍵參數

自激振蕩脈動特性主要由自激振蕩腔室結構和關鍵參數決定［7］，如圖2所示。關鍵結構參數如下：噴嘴進出口直管道管徑d；上噴嘴入口流道管徑d1；下噴嘴出口流道管徑d2；圓錐形收縮管夾角α1；下噴嘴碰撞壁夾角α2；圓錐形擴散管夾角α3；自激振蕩脈沖腔室長度L；自激振蕩脈沖腔室直徑DT。依據先前的實驗和仿真結果，自激振蕩腔室結構參數優選范圍和關鍵參數［13］，如表1所示。

圖2 自激振蕩腔室結構Fig.2 Self-Oscillating Chamber Structure

表1 自激振蕩腔室無量綱結構參數Tab.1 Self-Excited Oscillation Chamber Dimensionless Structural Parameters

3.2.2 模型驗證及網格無關性

計算方法選用大渦模擬，計算模型為可壓縮理想氣體，時間項為二階隱式格式，對流項和擴散項均為二階迎風格式。模型邊界條件，如圖3所示。根據先前實驗［16］，計算模型邊界條件采用壓力入口和壓力出口，上游入口壓力設為101325Pa，下游出口壓力設為84300Pa。

圖3 網格模型及邊界條件Fig.3 Grid Model and Boundary Conditions

由于3D-LES 計算成本較高，為了有效提升運算效率，采用2D 軸對稱大渦模擬是可行的［14-15］。通過與文獻［16］的實驗數據進行比較，測試了上游入口壓力為90500Pa條件下的平均流量，測試結果，如表2所示。采用2D-LES模型的誤差為1.41%，驗證結果表明數值計算結果和實驗結果之間具有良好的一致性。因此，該數值模型對于本研究是可靠的。

表2 2D-LES模型的平均流量誤差Tab.2 Average Flow Error of the 2D-LES Model

分別使用116282網格、194300網格及415721 網格，對自激振蕩腔室中心結構的2D-LES 計算模型進行網格無關性測試。軸向流速的網格無關性檢測，如圖4所示。綜合權衡計算效率和求解精度，最終確定選用194300網格模型。

圖4 網格無關性驗證Fig.4 Grid Independence Verification

4 多目標優化

通過面心復合設計、響應面法和NSGA-II算法進行多目標協同優化，獲取自激振蕩腔室無量綱結構參數的最優解，具體的多目標優化流程圖，如圖5所示。

圖5 自激振蕩腔室多目標優化流程圖Fig.5 Multi-Objective Optimization Flow Chart of Self-Excited Oscillation Chamber

4.1 約束處理

自激振蕩腔室減阻特性的關鍵要求是降低流阻，其中流阻主要包含兩部分，一部分是流體與固體壁面之間的摩擦阻力，另一部分是流體之間相對運動的內摩擦。因此，出口壁面平均摩擦阻力F1和出口流道平均內摩擦力F2作為多目標優化的兩個目標，計算公式如下：

式中：τ—單位面積上的剪切應力；A—自激振蕩腔室出口流道網格總面積；μ—內摩擦系數（黏滯系數）；du∕dy—速度梯度；V—自激振蕩腔室出口流道網格總體積。

4.2 試驗設計

中心復合設計（CCD）和面心復合設計（FCCD）是響應面分析的常用試驗設計方法之一。試驗設計的實驗點按照類別可以分為中心點、立方點及軸向點三種，不同類別實驗點位置示意圖，如圖6所示。

圖6 CCD和FCCD設計點分布對比Fig.6 Comparison of CCD and FCCD Design Points

可以看出在CCD 的設計點中，存在超出原定水平的數據。由于樣本點數據超出原定水平范圍違背了實際工況要求，因此，面心復合設計（FCCD）在本次優化中被采用，既保留了CCD試驗設計在擬合各個因素與輸出指標間的響應關系的優勢，又避免了樣本點數據超出原定水平。設計各因素水平及編碼值，如表3所示。

表3 試驗設計因素水平Tab.3 Trial Design Factor Level

4.3 響應面近似模型

響應面表面形式族的選擇對響應面的分析結果具有重大影響，因此，所選擇的響應面表面形式應滿足特定光滑度要求。對于不同的實際應用，最優的響應面表面形式的選擇也不相同。由于多項式的計算簡單，并且能夠生成封閉形式的表達式。在計算流體力學中，多項式成為表示響應面的常用形式。鑒于多目標優化的設計變量數為3，選擇二次多項式響應面的形式，其表達式為［12］：

式中：β0，βi，βij—回歸系數。響應面y表示為設計變量在設計空間內的函數。

在提出響應面近似模型用于多目標優化之前，應檢測響應面模型是否能夠滿足精度要求。如果不滿足精度要求，則需要增加樣本點，從而調整響應面的可變性和偏差之間的平衡，使得偏差減小，響應面可變性增強。R2準則和修正復相關系數R2adj通常被用于響應面的驗證，其表達式如下［12］：

式中：SSR—指回歸平方和，表示因回歸方程所引起的y的不均勻程度；SSY—指總體平方和，表示觀測值y的不均勻程度；SSE—誤差的平方和，表示因隨機誤差引起的y的不均勻程度；k—優化變量個數；n—樣本點個數。

4.4 多目標優化算法

在NSGA-II中，優化的最終結果是一組Pareto解集，也被稱為Pareto前沿。SBX作為交叉和突變的操作機制被用于生成下一代個體，包含交叉操作和變異操作，其表達式分別為［12］：

交叉操作：

變異操作：

式中：xi(1，t)和xi(2，t)—父代個體；xi(1，t+1)和xi(2，t+2)—子代個體；βqi—隨機變量；xiUB和xiLB—變量上下限；δq—變異算子；u—［0，1］之間的隨機數；ηm—指定的分布指數；δ—干擾因子。多目標優化算法的設置，如表4所示。

表4 NSGA-II優化算法參數設置Tab.4 NSGA-II Optimization Algorithm Parameter Settings

5 優化結果與分析

5.1 響應面模型構建

FCCD設置因素數為3，由于CFD模擬過程中不存在實驗環境等其他不可控因素影響，因此中心試驗重復次數為5，因變量個數為2，分別為出口壁面平均摩擦阻力F1和出口流道平均內摩擦力F2。為了消除兩個目標之間量綱的影響，對CFD 數值結果歸一化處理，試驗設計樣本點集及數值模擬所得結果，如表5所示。

表5 試驗設計和CFD歸一化結果Tab.5 Experiment Design and CFD Normalized Results

這里采用Quadratic 模型對目標變量F1，F2和設計變量d2∕d1，DT∕d1，L∕DT之間的函數關系進行擬合，將設計變量d2∕d1，DT∕d1，L∕DT分別記為x1、x2、x3，得到響應面回歸方程，如式（10）、式（11）所示：

5.2 模型顯著性分析及誤差檢驗

通過對響應面回歸模型進行檢測，方差分析結果，如表6所示。Prob＞F表示無顯著影響的概率，并且F值越大P值越小說明模型的相關性越顯著。當P小于0.01時，表示模型非常顯著。當P值大于0.01且小于0.05時，模型為顯著。當P值大于0.05表示模型不顯著。

表6 響應面回歸模型的方差分析結果Tab.6 Analysis of Variance by Response Surface Regression Model

對響應面回歸模型進行誤差統計分析，相關系數R2值和Radj2值越接近1，說明相關性越好。一般情況下，實際工程應用要求相關系數R2值和Radj2值大于0.9。響應面回歸模型誤差統計分析，如表7所示。統計分析結果表明，這里所得的響應面近似模型符合上述檢驗原則，具有較好的適應性。同時通過圖7可以看出，預測值與實際值基本分布在一條直線上，說明回歸模型具有較高的精度。

表7 響應面回歸模型誤差統計分析表7 Statistical Analysis of Errors in Response Surface Regression Models

圖7 預測值與實際值分布Fig.7 Distribution of Predicted and Actual Values

5.3 響應面關系及Pareto解集分析

響應面法能夠通過樣本點擬合出直觀的三維響應面，并且可以通過等高線圖的特征反映出設計變量之間的交互作用。通過控制某個因素不變的條件下，得到另外兩個因素的交互作用對出口壁面平均摩擦阻力F1和出口流道平均內摩擦力F2的影響。其中，等高線為橢圓形，表示交互作用顯著。d2∕d1，DT∕d1和L∕DT對出口壁面平均摩擦阻力F1的響應關系，如圖8 所示。可以看出，當L∕DT處于0 水平時，出口壁面平均摩擦阻力F1隨d2∕d1的增大而減小，隨著DT∕d1的增大而減小，并且無明顯交互作用。當DT∕d1處于0 水平時，出口壁面平均摩擦阻力F1隨著d2∕d1的增大而減小，并且隨L∕DT的增大出口壁面平均摩擦阻力F1先增大后減小，并且無明顯交互作用。當d2∕d1處于0 水平時，出口壁面平均摩擦阻力F1隨著DT∕d1的增大而減小，并且隨L∕DT的增大出口壁面平均摩擦阻力F1先增大后減小，并且無明顯交互作用。無論DT∕d1，L∕DT處在什么水平，F1隨著d2∕d1的增大都呈現出了明顯的降低趨勢，并且與DT∕d1，L∕DT相比d2∕d1起主導作用。從等高線圖可以看出并非橢圓形說明d2∕d1與DT∕d1，L∕DT對于出口壁面平均摩擦阻力F1的響應無明顯交互作用。d2∕d1，DT∕d1和L∕DT對出口流道平均內摩擦阻力F2的響應關系，如圖9所示。

圖8 d2∕d1，DT∕d1，L∕DT對出口壁面平均摩擦阻力的響應關系Fig.8 Response Relationship of d2∕d1，DT∕d1，L∕DTto the Average Frictional Resistance of the Outlet Wall

圖9 d2∕d1，DT∕d1，L∕DT對出口流道平均內摩擦阻力的響應關系Fig.9 Response Relationship of d2∕d1，DT∕d1，L∕DTto the Average Internal Frictional Resistance of the Outlet Flow Channel

由圖9可以看出，當L∕DT處于0水平時，出口流道平均內摩擦阻力F2隨d2∕d1的增大而先增加后減小，隨著DT∕d1的增大而增大，并且無明顯交互作用。當DT∕d1處于0水平時，出口流道平均內摩擦阻力F2隨著d2∕d1的增大而先增大后減小，隨L∕DT的增大而先增大后減小，并且從等高線投影結果可以看出有明顯交互作用，出口流道平均內摩擦阻力F2存在極大值點。當d2∕d1處于0水平時，出口流道平均內摩擦阻力F2隨著DT∕d1的增大而增大，并且隨L∕DT的增大出口流道平均內摩擦阻力F2先增大后減小，并且無明顯交互作用。

NSGA-II 算法結合響應面模型進行多目標優化計算，得到Pareto最優前沿，如圖10所示。由于在Pareto前沿中，沒有任何一個解能夠優于其他解，因此任何一個解都是可以被選擇Pareto最優解。

圖10 目標F1和F2的Pareto最優前沿Fig.10 Pareto Front for Objectives F1 and F2

5.4 中心結構與優化結構數值模擬對比

通過多目標優化，充分考慮加工精度和設計標準，引入1：1的決策權重，從Pareto前沿中選出1個符合設計條件的候選設計點進行CFD仿真，并與中心設計點原始結構參數進行對比，結果如表8所示。

表8 中心結構與優化結構對比Tab.8 Comparison of Central Structure and Optimized Structure

從優化結果可以得出，優化后的模型與CFD 模擬值之間存在一定的誤差，分別為2.18% 和4.58%，滿足精度要求。與中心原始結構相比，出口壁面平均摩擦阻力和出口流道平均內摩擦力分別減少了43.01%和54.52%。

6 結論

這里對自激振蕩腔室的無量綱結構參數進行了多目標優化，以獲得最優減阻效果。采用了2D大渦模擬的數值分析方法，結合面心復合設計（FCCD）、響應面法和NSGA-II算法進行了多目標協同優化。分析了不同無量綱結構參數對于自激振蕩腔室減阻特性的影響。得到如下結論：

（1）通過CFD與面心復合設計（FCCD）、響應面法和NSGA-II算法的多目標協同優化方法，能夠有效獲取自激振蕩腔室無量綱結構參數的最優解，實現自激振蕩腔室的最佳減阻。

（2）分析響應面關系，得出了d2∕d1，DT∕d1，L∕DT與出口壁面平均摩擦阻力F1和出口流道平均內摩擦力F2之間的關系。出口壁面平均摩擦阻力F1隨著d2∕d1的增大都呈現出了顯著降低，并且與DT∕d1，L∕DT相比d2∕d1起主導作用，d2∕d1，DT∕d1，L∕DT無明顯交互作用。出口流道平均內摩擦阻力F2s隨著DT∕d1的增大都呈現出了顯著升高，并且與d2∕d1，L∕DT相比DT∕d1起主導作用，d2∕d1與L∕DT之間存在明顯交互作用，并存在F2極大值點。

（3）充分考慮加工精度和設計標準，引入1：1的決策權重，從Pareto前沿中選出1個符合設計條件的候選設計點作為優化后的自激振蕩腔室無量綱結構。與中心原始結構進行對比，結果表明，優化后的自激振蕩腔室減阻性能得到了顯著提高。與中心原始結構相比，出口壁面平均摩擦阻力和出口流道平均內摩擦力分別減少了43.01%和54.52%。