陳蔣力 陳少強 任峰? 胡海豹?
1) (西北工業大學航海學院,西安 710072)
2) (中國船舶重工集團公司第七〇五研究所,西安 710077)
針對經典圓柱繞流問題,采用深度強化學習方法,提出了基于壁面壓力反饋的圓柱繞流減阻閉環主動控制方法,并比較分析了施加控制前后圓柱阻力系數、升力系數及流場的差異.控制系統中,以圓柱壁面上均勻分布的壓力探針測得的信號作為反饋,利用多層感知機建立壓強信號與吹/吸射流及控制效果的映射關系,即控制策略;通過在圓柱上下表面狹縫施加連續可調的吹/吸射流來進行主動控制.同時,利用深度強化學習中的近端策略優化方法,在大量的學習過程中對該控制策略進行不斷調整和優化,以實現穩定減阻效果.在圓柱繞流流動環境搭建方面,采用格子Boltzmann 方法建立與深度強化學習模型之間的交互式框架,模擬提取非定常流場條件下圓柱表面的壓強信號,并計算實時調整吹/吸射流強度時圓柱表面升力、阻力數據,以評估所選控制策略的優劣.研究表明:雷諾數為100 時,主動控制策略能減少約4.2%的圓柱阻力,同時減少約49%升力幅度;同時施加主動控制后圓柱的減阻效果與圓柱回流區長度呈現強相關趨勢.此外,不同雷諾數下智能體習得的策略減阻效果不同,雷諾數為200 和400 時,該主動控制策略能依次減小圓柱阻力17.3%和31.6%.本研究可為后續開展基于壁面壓力反饋的圓柱流動主動控制實驗以及其他復雜環境下鈍體流動智能控制提供參考.
流動控制根據是否需要外部能量輸入分為被動流動控制和主動流動控制[1].被動控制具有無需能量輸入、設置簡單、成本低等優點,但當流場實際情況偏離預期時,控制往往難以達到最佳效果[2].在主動控制中,根據是否需要從流場中獲取反饋信息并用以調節激勵器的輸出大小,可以分為開環控制和閉環控制[3].研究發現,閉環主動控制相較于開環控制具有自適應能力強的特點,能夠充分發揮激勵器的潛力,在少量能量輸入情況下往往就能得到很好的控制效果,因而有著更大的潛力[4].然而,非線性的Navier-Stokes 方程所描述的流場具有高維、多模態等復雜特征,導致設計有效的實時閉環主動流動控制策略往往十分困難[5].
近年來,迅速發展的深度強化學習(DRL)以其擅長與環境進行交互的特點,為上述流動控制難題帶來了新的思路.研究表明,深度強化學習能夠在高維、非線性等復雜環境中有效地習得控制策略[6].如果利用深度強化學習與流動控制環境進行交互,在不斷試錯、調整優化策略后,由習得策略建立控制律,對于閉環流動控制方法來說具有重要意義.Rabault 等[7]將深度強化學習應用于鈍體減阻,利用近端策略優化方法[8]成功習得了能夠實現穩定減阻約8%的閉環主動控制策略.該研究中,將位于圓柱周圍和下游流場中的151 個傳感器(每個傳感器同時采集流向速度及橫向速度)測得的速度作為反饋信號,并將獎勵函數(與升力系數和阻力系數有關)作為獎懲機制.Paris 等[9]引入一種新算法(S-PPO-CMA)優化傳感器位置,在雷諾數120 時,成功習得了穩定減阻約18.4%的閉環主動控制策略.任峰等[10]基于格子Boltzmann 方法搭建流動環境,在弱湍流條件下成功習得減阻的閉環主動控制策略,實現了約30%的圓柱減阻效果,為優化傳感器布局,進行了敏感性分析.
上述研究均以尾流中的速度傳感器等作為反饋信號.考慮到應用條件下的實際情形,在尾流固定位置處布置傳感器技術實現上比較困難,因此本研究以實際應用場景中較容易測得的壁面壓力作為反饋信號.本文采用GPU 加速的格子Boltzmann方法[11]對低雷諾數下的圓柱繞流進行數值模擬,將深度強化學習應用于主動流動控制,并通過分析施加控制前后射流速度的變化得到主動控制對流動的影響機制.
本文選取了圖1 所示的物理模型.其中,直徑為D的圓柱放置在長L=21.79D,寬H=4.06D的狹窄管道中.將圓柱中心設置為坐標系零點,與入口邊界距離L1=2D,與出口邊界距離L2=19.79D.為了促進渦旋脫落,圓柱中心略微偏離計算域中心線(y方向上向上偏離0.05D).由于計算域較為狹窄,上下兩側壁面會對圓柱的升阻力系數有所影響,但不影響最終控制效果.雷諾數定義為Re=UmD/υ,其中D為圓柱直徑,υ是流體的運動黏度,Um為入口處的平均速度.參考時間T定義為T=D/Um.

圖1 物理模型示意圖Fig.1.Schematics of the physical model.
本文圓柱的阻力系數和升力系數分別用CD和CL表示,定義如下:


其中ρ為流場中流體的密度,FD以及FL分別為圓柱所受的總阻力和總升力.
本文通過圓柱上下表面的狹縫施加吹/吸射流控制,狹縫寬度與圓心角為10°的弧長相對應,為保證零質量流量射流,上下表面射流速度的大小與方向均相同.為了保證控制的連續性,并避免施加控制后射流速度的跳躍,將射流速度表示為

其中Ujet與分別為當前時間步以及下個時間步的射流速度,兩個時間步的時間間隔為T/δt=800 (δt=1 為格子時間),at指在當前時間步下,深度強化學習智能體輸出的動作值(智能體及動作值的詳細解釋在2.4 節給出).α為數值參數,α=0.1 與Rabault 等[12]一致.(3)式為指數衰減形,其表達形式及α的取值會影響訓練效率,詳見參考文獻[12].
本文獎勵函數設置為

其中〈·〉 代表在兩次施加激勵過程中進行時間平均,w為權重,本文將權重設置為1,這樣不僅能減少圓柱阻力,還能保證在較低雷諾數時圓柱升力系數不會出現顯著零偏,詳見參考文獻[10].
由于格子Boltzmann 方法不依賴連續介質假設,算法簡單,具有良好的計算局部性等優點[13,14].本文基于格子Boltzmann 方法建立流動數值模擬環境.本文的格子Boltzmann 求解器使用了D2Q9格式[15](二維流動,每個格點處速度向9 個方向離散)作為離散速度模型,通過多松弛時間算法[16,17]來提高數值模擬的精度和穩定性,采用He-Luo 模型[18]確保本文中流體的不可壓縮性.
本文中多松弛時間模型表達式為

式中f(x,t)指在t時刻x位置處的分布函數,feq為平衡態分布函數,ei為粒子的離散速度,δt為格子時間步長,M稱為變換矩陣,可以將速度空間轉換至矩空間.松弛參數設置為

其中τ為松弛時間,松弛時間與運動黏度關系為τ=0.5 +,cs=3–1/2為格子聲速.
邊界條件的設置與此前工作[10]類似,入口處的速度設置為拋物線形,數學表達式為

出口處的壓力設置為零,入口和出口均采用非平衡外推格式[19].上下壁面采用具有二階精度的半步長反彈格式[20],圓柱曲面使用雙線性插值法進行處理[21],上下壁面及圓柱表面均為無滑移壁面.此外,圓柱受力采用改進的動量交換方法求得[22].
網格數量以及時間分辨率的選取,對數值模擬時的計算精度和計算耗時有一定的影響,表1 測試了Re=100 時三種不同網格數和時間分辨率下的算例,對比了這些算例中圓柱的各個參數.其中T為參考時間,Sr為斯特勞哈爾數(Sr=f D/Um,其中f為渦旋脫落頻率).

表1 無關性驗證(Re=100)Table 1.Validation and convergence study (Re=100).
從表1 中可以看到,與算例Ⅲ(本文最佳的網格數及時間分辨率)相比,算例Ⅰ中平均阻力系數與平均升力系數的絕對值分別相差0.009 和0.028,相對誤差分別為0.3%和4.4%.而算例Ⅱ的平均升力系數的絕對值相對誤差在0.5%以內,且計算耗時小于算例Ⅲ,很好地滿足了低計算耗時和高計算精度的要求.因此,Re=100 時,本文以表1 中的算例Ⅱ為標準,搭建流動數值模擬環境.
在合適的時間分辨率和網格數下,成功搭建了Re=100 時的流動數值模擬環境.表2 給出了該環境下的計算結果以及其他學者數值模擬的結果,包括圓柱表面的最大阻力系數、最大升力系數以及渦旋脫落頻率(Sr).對比發現,各個計算結果相對偏差均在2.5%以內,這驗證了格子Boltzmann 方法在該問題上的準確性與可靠性.

表2 Re=100 時最大升阻力系數CD,max,CL,max 與Sr 對比Table 2.Comparison of CD,max,CL,max and Sr at Re=100.
深度強化學習的本質是互動學習[25].智能體(本文指模擬人的思維決策過程,可以同環境相交互的程序)主要通過以下三個方面與環境進行交互:動作at,狀態st和獎勵函數rt.本文所使用的是深度強化學習中的近端策略優化方法(PPO),具有較高的穩定性和可靠性等優點[8].
首先引入軌跡τt的概念,表示在給定參數θ下,某次學習時,狀態-動作-獎勵函數隨時間的變化關系:

(8)式指環境輸出狀態值s0給DRL 智能體,選擇動作a0,并得到獎勵函數r0,環境被動作改變,輸出新的狀態s1.經過不斷重復,得到某次學習下的控制策略,并以累計獎勵函數R(t)評估策略,其中γ(0 <γ<1)為折扣系數,本研究中固定為0.97.
近端策略優化方法使用了兩組人工神經網絡,分別為動作策略(輸入狀態,輸出動作的概率分布)的actor 網絡和評價策略(輸入狀態,得到值函數,計算出優勢函數)的critic 網絡,即actor-critic網絡.創建critic 網絡用以減小實際獎勵與預期獎勵的數學期望之間的差異,即


其中參數θ代表神經網絡的所有權重的集合,V(st)指值函數,即狀態st下所有可能的動作與其累計獎勵函數的乘積的和(預期獎勵的數學期望).
本文依據Schulman 等[8]的工作,選取了近端優化策略中actor 網絡的目標函數,用于輸出某狀態下動作的概率分布.該目標函數算法簡單且具有良好的穩定性[8].其數學表達形式為:

其中,rt(θ)指的是新舊策略概率比(πθ(at|st) 和πθold(at|st)分別為當前策略和舊策略在狀態st下采取動作at的概率).clip 函數用于限制新舊策略概率比,是指將rt(θ)的大小限制在[1–ε,1+ε]之間,當rt(θ)的值大于1+ε時,用1+ε代替rt(θ),當rt(θ)小于1–ε,用1–ε代替rt(θ).ε是超參數,設置為0.2,與Schulman 等一致,在Schulman 文中,與ε=0.1 和ε=0.3 相比,ε=0.2 時算法的效果最佳,詳見參考文獻[8].min()即取下界,具體理解為:當>0 時,當前動作產生的獎勵估計大于預期獎勵的數學期望,因此增加新策略中當前動作出現的概率,為提高訓練的穩定性,限制其不能超過原策略的1+ε倍;反之,當<0 時,限制其不能小于原策略的1–ε倍.
本節主要復現了Rabault 等[7]的結果,并采用了相同位置的速度傳感器.Re=100 時,智能體以流場速度作為反饋信號,訓練得到減阻的閉環主動控制策略.由于已有學者[7,10,12]開展過相似的研究內容,此處不再詳細描述.
圖2 給出了基于速度反饋的智能體的訓練過程,本文進行三次不同的訓練,并對三次訓練過程中的阻力系數取平均值,得到一條更為平滑的曲線.每一訓練集(episode)均相當于一次完整的數值模擬過程,持續時間為32T,可以看到,隨著訓練集數的增加,在200 集后,阻力系數(通過對每一集后半段的阻力系數取平均得到對應阻力系數)逐漸穩定,在400 集時阻力系數大約在2.98 附近,但仍然有大幅度波動出現,這是由于深度強化學習智能體的訓練過程具有一定的隨機性.

圖2 智能體訓練過程(速度反饋)Fig.2.Learning curves of DRL agent (velocity feedback).
圖3 為訓練完成,智能體已習得閉環主動控制策略后,施加控制前后圓柱的阻力系數隨時間變化對比圖,可以看到,施加控制至流場達到穩定后,圓柱的阻力系數為2.989,減阻率為6.6%,且射流速度僅為入口處平均來流速度的0.15 倍.

圖3 控制前后阻力系數(CD)隨時間變化圖Fig.3.Temporal variations of drag coefficient (CD) with and without active flow control.
在成功搭建流動控制環境后,智能體需要基于有限的流場信息來訓練主動流動控制策略.考慮到速度傳感器成本大,在尾流處布置困難,難以用于實際實驗中.而壓力傳感器較速度傳感器成本低,在圓柱表面處布置容易.本節通過壁面壓力反饋方法代替原來的尾流速度反饋方法.
考慮到壓力探針數量對智能體的影響,本文分別對比了6,14 和30 個壓力探針下智能體習得控制策略的減阻效果.圖4 分別給出了6,14 及30 個壓力探針的位置分布,采用在壁面上均勻分布的方式并避開了圓柱上下表面的吹/吸射流.

圖4 壓力探針位置圖 (a) 6 個壓力探針;(b) 14 個壓力探針;(c) 30 個壓力探針Fig.4.Schematics of the pressure sensors position:(a) 6 pressure sensors;(b) 14 pressure sensors;(c) 30 pressure sensors.
本文圓柱曲面邊界處理方法為雙線性插值法[21],由于圓柱表面位于流體格點與固體格點之間,無法直接得到壓力分布曲線.因此采用了多次線性插值的方法,通過圓柱表面附近格點間接得到未控制下圓柱的壓力分布曲線.
為了和Tiwari 等[26]進行對比,本節搭建了Re=20 時的數值模擬環境,并得到圓柱表面的壓力分布曲線.其中,壓力經過無量綱處理,即

式中,Pref為參考壓力,這里Pref=1907,Umax=1.5Um是入口處的最大速度.
Re=20 時圓柱表面壓力分布曲線對比圖如圖5 所示,其中點(0.5,0)作為橫坐標起點,逆時針方向為正方向.可以看到本文得到的Re=20時的圓柱壓力分布曲線與Tiwari 等[26]的結果幾乎重合.

圖5 圓柱壓力分布曲線,Re=20Fig.5.Pressure distribution curve along the cylinder surface,Re=20.
Protas 和Wesfreid[27]將圓柱的阻力分為兩部分,一部分為由渦旋脫落引起的阻力,另一部分為圓柱無渦旋脫落時的阻力.研究表明,對于圓柱繞流而言,主動控制僅能減少由渦旋脫落引起的阻力[27?29].同時,圓柱阻力與回流區長度具有一定的聯系,在無渦旋脫落時圓柱尾流處回流區呈細長狀,長度隨著Re的提高而增加.施加主動控制后,圓柱阻力減小,此時回流區長度更接近穩定無渦旋脫落時的回流區長度.Rabault 等[7]指出,Re=100 時,本文計算域下無渦旋脫落時圓柱的阻力系數約為2.93.
本節智能體以壁面壓力作為反饋信號,并進行訓練,訓練曲線如圖6 所示.與圖2 相比,圖6 中的智能體需要更多集數進行訓練以使阻力系數趨于穩定.訓練完成后,將訓練過程中阻力系數最小時對應的控制策略作為本次訓練得到的最佳策略.智能體分別采用6,14 和30 個壓力探針作為反饋信號,通過主動控制和無控制得到的圓柱阻力系數、升力系數和射流速度隨時間變化曲線,如圖7所示.

圖6 智能體訓練過程(壓力反饋)Fig.6.Learning curves of DRL agent (wall pressure feedback).

圖7 主動控制和無控制下圓柱的阻力系數(CD),升力系數(CL)和射流速度(Ujet)變化曲線,Re=100Fig.7.Temporal variations of drag coefficient (CD),lift coefficient (CL),and jet velocity (Ujet) with and without active control,Re=100.
如圖7 所示,探針數量不同的智能體施加主動控制均能減少圓柱阻力.當壓力探針數量為6 時,能習得減阻的閉環主動控制策略,但是僅將圓柱的平均阻力系數減少到3.11 左右,平均升力系數由–0.026 增加至–0.019,減阻效果達到了2.9%,減少的阻力占因渦旋脫落產生的阻力的34.1%.相比之下,壓力探針數量增加至14 時,智能體習得的控制策略能將平均阻力系數減小到3.07 左右,平均升力系數由–0.002 減少至–0.013,減阻效果達到了4.2%,減小的阻力占由于渦旋脫落導致的阻力的51.3%.當壓力探針增加至30 個時,習得的控制策略也只能將平均阻力系數減小到3.08 左右,平均升力系數由0.010 增加至0.012.同時,升阻力系數的變化幅度減少.阻力波動從無控制下的[3.16,3.25]分別減少至[3.08,3.13](對應6 個壓力探針),[3.05,3.09](對應14 個壓力探針),[3.06,3.09](對應30 個壓力探針).升力波動從無控制下的[–0.99,1.01]分別減少至[–0.62,0.63](對應6 個壓力探針),[–0.53,0.52](對應14 個壓力探針),[–0.55,0.57](對應30 個壓力探針).從射流速度隨時間變化曲線中可以看出,主動控制時,首先以較大的射流速度改變流場結構,此時與圓柱阻力快速減少相對應.在控制達到穩定后,僅需少量能量(以14 個探針為例,約為剛開始控制時最大射流速度的0.5 倍,入口處平均速度的0.15 倍)即可以達到很好的減阻效果.
本節試圖分析控制前后流動的變化,來探究主動控制對流場的影響.圖8 為控制前后的瞬時流場云圖(以14 個壓力探針習得的控制策略為例),每張圖的上部分為無控制下的流場云圖,下部分為施加主動控制下的流場云圖.由于在Re=100 時進行數值模擬,云圖中可以看到規則的流場,受逆壓梯度和流體黏性的影響,圓柱壁面產生邊界層分離現象,并在分離點后形成回流區(流向速度為負的區域)和尾流區(圓柱尾渦周期性脫落).

圖8 Re=100 時的瞬時流場云圖(a1)—(d1)無控制下流向速度,橫向速度,壓力及渦量云圖;(a2)—(d2)主動控制下流向速度,橫向速度,壓力及渦量云圖Fig.8.Instantaneous contours of flow fields at Re=100:(a1)–(d1) Contours of streamwise velocity,transverse velocity,pressure,and vorticity without active control;(a2)–(d2) contours of streamwise velocity,transverse velocity,pressure,and vorticity with active control.
從瞬時的流向速度對比云圖中可以看到,與控制前相比,施加主動控制后圓柱上分離點位置變化較小,回流區長度顯著增加,達到了2.08D左右,而控制前回流區長度僅有1.68D左右,回流區長度增加了23.8%.
從瞬時橫向速度對比云圖和瞬時壓力對比云圖中可以看到,相比無控制下的流場,主動控制下流場內流體的瞬時橫向速度變小,圓柱前后的壓力差(圓柱的壓差阻力)減小,尾流處壓力為負的區域面積增大.從瞬時渦量對比云圖中可以看到,主動控制下尾流處渦的形成和脫落過程被改變,渦開始脫落的位置和未控制時相比向下游推移,同時渦脫落頻率和渦量強度減小.
對500 張云圖(約3 個渦脫落周期)取平均得到如圖9 所示的時間平均的流向速度云圖、橫向速度云圖和壓力云圖.分析流場時均云圖同樣看到,與無控制情況相比,施加主動控制后圓柱分離點位置變化較小,回流區長度顯著增加,時均壓力云圖中,圓柱前后的壓力差變小,尾流處壓力為負的區域面積增大,這說明施加主動控制能減少圓柱的部分壓差阻力.

圖9 Re=100 時的時均流場云圖 (a1)—(c1)無控制下流向速度,橫向速度及壓力云圖;(a2)—(c2)主動控制下流向速度,橫向速度及壓力云圖Fig.9.Time-averaged contours of flow fields at Re=100:(a1)–(c1) Contours of streamwise velocity,spanwise velocity,pressure,and vorticity without active control;(a2)–(c2) contours of streamwise velocity,spanwise velocity,and pressure with active control.
在第2 節中,智能體始終在Re=100 時的數值模擬環境中進行訓練.隨著Re的改變,繞流流場也會發生變化,研究指出當Re增加時,因渦旋脫落產生的阻力占總阻力的比例提高,且無渦旋脫落時圓柱尾流處回流區長度也隨之增加[25,27?29].本節通過對比不同Re智能體習得的閉環主動控制策略(以第3 節14 個壓力探針為例),探究Re對減阻效果及流場的影響.
Re=100,200 和400 時,智能體施加主動控制后得到的結果如圖10 所示.圖10 中的數據表明,Re不同時,智能體均能習得減阻的閉環主動控制策略,在本文研究范圍內(Re=100—400),無控制下圓柱的阻力系數與升力系數的變化幅度隨著Re的增大逐漸提高,主動控制下圓柱減阻率同樣隨著Re增大逐漸提高.當Re=100 時,主動控制下圓柱的平均阻力系數達到穩定時僅僅只有3.07,阻力系數變化幅度從未控制下的0.096 減少至主動控制下的0.054,減阻效果達到了4.2%;當Re=200 時,主動控制下圓柱的平均阻力系數達到穩定時為2.59,減阻效果達到了17.3%;Re=400 時,施加控制后圓柱的平均阻力系數最終穩定在2.15 左右,阻力系數幅值從未控制下的0.504 減少至主動控制下的0.129,減阻效果高達31.6%.不同Re下升力系數變化幅度均減少,Re=100 時,升力系數變化幅度由[–0.98,1.01]減少至[–0.53,0.52];Re=200 時,升力系數變化幅度[–2.11,2.16]由減少至[–1.09,1.22],控制下升力系數平均值為–0.025;Re=400 時,升力系數變化幅度[–2.78,2.86]由減少至[–0.73,0.33],控制下升力系數平均值為–0.212.

圖10 主動控制與無控制下的圓柱阻力系數(CD),升力系數(CL)和射流速度(Ujet)變化曲線Fig.10.Time-resolved value of drag coefficient (CD),lift coefficient (CL),and jet velocity (Ujet) with and without control.
通過研究圖10 中射流速度隨時間變化的曲線,可以發現閉環主動控制策略的具體實施方案.如第3 節所述,施加主動控制時,智能體首先以較大的射流速度改變流場結構,隨后,以較小的射流速度維持流場形態.和Re=100 時的射流速度不同,Re=400 時射流速度在達到平衡時出現顯著的零偏,本文認為這是因為當Re增加后,圓柱減少了更多的阻力,使得阻力系數的絕對值在獎勵函數中的占比較大,這使得升力系數可能出現較大的零偏.可以通過提高獎勵函數中升力系數的權重以減少射流速度的零偏,詳見參考文獻[10].
Re=200 和Re=400 時控制前后的時均流向速度云圖如圖11 所示.可以發現,在本文的Re范圍內,圓柱減阻率與圓柱回流區長度強相關,Re越高,主動控制下圓柱減阻率越大,回流區長度越長.Re=200 時,回流區長度增加136.7%;Re=400 時,回流區長度增加341.5%.

圖11 無控制下和施加主動控制下的時均流向速度云圖(a1) Re=200,無控制;(a2) Re=200,施加控制;(b1) Re=400,無控制;(B2) Re=400,施加控制Fig.11.Time-averaged streamwise velocity fields without control and with active control:(a1) Re=200,without control;(a2) Re=200,with control;(b1) Re=400,without control;(b2) Re=400,with control.
開展了基于壁面壓力的圓柱減阻智能流動控制研究,得到以下結論.
1) 在智能體進行學習的過程中,壓力探針數量的選取較為重要.壓力探針數量太少,會影響減阻效果,壓力探針數量太多,則會增加控制系統成本.
2) 主動控制時,圓柱的阻力減少,對應圓柱的回流區長度顯著增加,渦強度減小.同時Re=100 時圓柱的升力系數波動幅度減少49%.
3) 主動控制時,先以較大的射流速度改變流場結構,達到穩定后,以較小的射流速度維持流場結構,即只需較小的能量輸入可實現較好的減阻效果.
4)Re不同時,智能體習得的閉環主動控制策略不同,在本文所選Re范圍內,主動控制時圓柱的減阻率隨Re增加而增加,Re=100 時,圓柱的減阻率達到了4.2%.Re=200 時,圓柱的減阻率達到了17.3%.Re=400 時圓柱的減阻率達到了31.6%.