基于流場反演和圖神經網絡的翼型分離流動預測方法

摘要： 傳統的RANS模型采用布辛涅斯克近似，假設湍流雷諾應力和平均速度梯度張量之間呈線性關系，這一假設適用于簡單的剪切流動，但很難推廣應用于復雜分離流動問題. 本文基于流場反演和機器學習FIML 方法框架，針對目前該方法框架內普遍采用的多層感知機網絡對于湍流空間相關性表征不足的缺陷，通過圖神經網絡的引入，對工程應用較為廣泛的SA一方程湍流模型生成項進行了修正. 在此基礎上，結合流場分離特征設計了一種加權函數，改進了圖神經網絡的消息傳遞機制. 針對大攻角、高雷諾數S809 翼型分離流動的實驗結果表明，與現有的多層感知機網絡相比，圖神經網絡在不同攻角、不同網格上預測得到的升力系數和實驗值更加接近，且新的消息傳遞機制能夠進一步提升圖神經網絡預測精度.

關鍵詞： 湍流建模； 流場反演； 圖神經網絡； 計算流體力學； 數值模擬

中圖分類號： TP391 文獻標志碼： A DOI： 10.19907/j. 0490-6756. 240008

1引言

湍流的精確模擬對于認識復雜流動機理和實現復雜工程應用至關重要［1］. 常用的湍流模擬方法包括直接數值模擬、大渦模擬和雷諾平均模擬等.直接數值模擬（Direct Numerical Simulation，DNS）和大渦模擬（Large Eddy Simulation， LES）能夠對一些簡單流動進行精準的預測，但是對于高雷諾數下的復雜流動，計算開銷和時間成本極高，當前的高性能計算資源很難支撐此類方法在工程外形上的應用［2，3］. 目前工程湍流問題的模擬主要依賴雷諾平均模型（ 例如k- ε，k- ω，Spalart-Allmaras（ SA），Shear Stress Transport（ SST） 和雷諾應力模型）. 但是，RANS（Reynolds-averagedNavier-Stokes）模型采用布辛涅斯克假設（渦粘性假設），通過求解湍流模型方程得到渦粘性系數，進而完成對于流動控制方程的封閉. 但RANS 模型中的參數通常采用簡單的流動問題進行標定，這也使得RANS 模型在預測逆壓梯度、分離等復雜流動現象時會產生較大的誤差［4］. 因此，為了提高湍流模擬的精度，構建一個能夠適應多種流動狀態的湍流模型具有重要意義.

近年來，利用機器學習算法（Machine Learning，ML）和高保真的實驗數據，數據驅動的方法越來越多地應用于傳統計算流體力學（ComputationalFluid Dynamics， CFD）領域. 從圖網絡（Graph Network）［5］的提出，到Sanchez-Gonzalez等［6］首次將圖網絡應用于物理問題的求解與模擬，再到Sanchez-Gonzalez 等［7］將圖網絡結合粒子基離散化方法來模擬水的運動，圖網絡模型在物理模擬研究領域表現出強大的學習能力. 湍流建模作為CFD 領域中的一個重要課題，數據驅動湍流建模的研究主要分為以下兩類方法.

（1）傳統湍流模型修正. 這類方法認為造成流場計算的偏差主要源于湍流模型方程的經驗系數［8］. Cheung 等［9］應用貝葉斯推理對SA 湍流模型的不確定性進行量化，并對參數進行校正. 同樣，Ray 等［10，11］研究了k- ε 湍流模型中的參數不確定性，并將校準參數應用于交叉流中的射流. Duraisamy等［12］總結了校準模型參數的進一步研究.不足的是，對于不同類型的流動，校準后的模型參數并不完全相同. 為了克服這一缺陷，Duraisamy等［12］和Singh 等［1，13，14］提出了一種數據驅動建模框架，稱為流場反演和機器學習（Field Inversion andMachine Learning，FIML）. 該方法通過基于伴隨的方法獲得空間修正系數對RANS 模型進行了修正. 然后，利用機器學習算法來預測修正系數，并在各種相似的流動中獲得了改進的預測效果.Yang 等［15］采用正則化集成卡爾曼濾波的方法來解決基于流場反演和機器學習FIML 框架中的反問題，降低了代碼修改的復雜性. 不同的是，在上述的研究中機器學習訓練采用離線訓練模式，Holland等［16］將機器學習算法集成到了流場反演過程中，將兩者耦合起來，采用在線訓練模式，提高了模型的泛化能力.

（2）完全數據驅動模型［4］. 這種方法不是基于現有的湍流模型，而是根據求解器求解出的流場變量的物理意義，構造具有一定物理意義的特征作為機器學習的輸入和輸出，利用機器學習算法對輸入特征與輸出特征之間潛在的關系進行建模，建立輸入特征與輸出特征之間的映射關系，由此便得到了新的模型，輸入特征與輸出特征之間的關系是高度非線性的. Zhu 等［17，18］和Sun 等［19］建立了由SA 模型數據訓練的數據驅動的線性渦粘模型，該方法在數據驅動模型和RANS 求解器之間采用了迭代模式. 后續的研究表明，相比于凍結耦合模式，迭代模式能夠有效地提高收斂性，特別是在高雷諾數情況下效果更佳. 相比于線性渦粘模型，非線性渦粘模型的適用范圍更加廣泛，數據驅動非線性渦粘模型的一個典型框架是張量基神經網絡（Tensor Basis Neural Network，TBNN）. 受Pope［20］工作的啟發，即將雷諾應力表示成了10個基張量的線性組合，Ling 等［21，22］提出了TBNN，利用深度神經網絡，針對雷諾應力的各向異性建立張量基神經網絡模型，并對雷諾應力分量進行計算，并將預測得到的雷諾應力代入RANS 求解器，得到了改進的速度場. Ling 等［23］根據流場變量的物理意義，提出了12 個具有旋轉不變性的變量作為機器學習的輸入特征，使得結果得到極大的提高，因此得出結論：輸入特征的旋轉不變性是模型輸入特征的選取原則. Wu 等［24］將流場的應變率張量和渦量張量梯度等進行組合，形成了較為完備的輸入特征組，使得模型預測能力進一步提升，這也表明輸入特征的完備性也是輸入特征的選取原則之一. Zhang 等［25］采用集成卡爾曼濾波法，通過間接觀測數據對TBNN 進行訓練. 除此之外，Weatheritt 等［26，27］利用基因表達編程建立了非線性渦粘模型的顯式表達.

總的來說，對于第一類方法，即對傳統湍流模型的修正方法，這類方法通常采用基于伴隨的方法和集成卡爾曼濾波法進行流場反演，從而得到修正系數的空間分布，隨后利用隨機森林（RandomForest）和深度神經網絡（Deep Neural Networks，DNN）等機器學習方法構造流場變量與修正系數之間的映射函數. 相比于第二種完全數據驅動模型的方法，傳統湍流模型已經被廣泛使用數十年，在多種流動情況下的效果也已經得到驗證，并且傳統模型通常基于流體力學基本原理，自然包含了一些物理約束和先驗知識，確保模型預測在物理上是可行的. 因此第一類方法目前被學者們廣泛研究，Wu 等［24］也指出目前機器學習的主要作用應是對現有湍流模型進行改進，以降低流場計算時的誤差，而不是徹底的取代湍流模型.

因此，本文基于第一類方法，將FIML 方法應用于本文的研究當中. 同時，在傳統FIML 框架中大多采用MLP 進行建模，在MLP 這種點對點的映射中，并沒有考慮湍流這種真實世界復雜物理現象的空間關聯性. 針對目前該方法對于湍流空間相關性表征不足的缺陷，介于GNN 所具備的空間關系建模能力，本文提出提出使用圖神經網絡進行湍流建模. 搭建了一種基于Encoder-Process-Decoder 結構的圖神經網絡模型，引入了網格點的相對位移向量以及兩點間的歐氏距離等空間信息，使得多個點能夠同時參與預測一個點上的SA湍流模型修正系數，并且提出了一種新的訓練策略，證明了在FIML 框架中引入空間信息的有效性，在高雷諾數的S809 翼型算例上得到了驗證.

2本文方法及模型

針對目前SA 湍流模型對翼型分離流動預測不準確的問題，本文基于流場反演和機器學習（FIML）框架，在SA 模型中引入一個修正項進行改進. 考慮目前CFD 求解所用到的網格正好適用于圖這種數據結構，提出將流場反演得到的結果與圖神經網絡結合，并通過觀察求解得到的特征分布，設計了一種基于特征分布的訓練策略，以高精度的數據驅動訓練，從流場中提取物理參數訓練模型，達到對湍流模型的優化，與當前大多使用的多層感知機網絡（MLP）相比取得了更好的效果.

2. 1流場反演與機器學習框架

將FIML 方法應用于SA湍流模型，以提升其準確預測翼型分離流動的能力. SA 湍流模型的控制方程的簡潔表達形式如式（1）

其中，P，D，T 分別表示生成項、破壞項以及擴散項. 由于SA 湍流模型計算復雜流動并不準確，因此引入修正項后表達式如式（2）.

引入修正項β（x）之后，不僅是修改了生成項，并且影響了整個源項的平衡. 接下來，根據FIML方法來預測修正項β（x），改進SA 模型對翼型分離流動的預測效果. FIML 方法的整體框架如圖1 所示，主要由3個關鍵步驟構成：（1） 流場反演，用實驗所獲的高保真參考數據去校正RANS 模型的求解結果，不斷優化得到修正項β（x）的空間分布.（2） 模型訓練，利用機器學習算法，訓練圖神經網絡模型，模型輸入是從流場中提取到的包含物理信息的流場特征，輸出是通過流場反演得到的修正項，構建從流場特征到修正項的映射.（3） 模型預測，利用步驟（2）訓練的模型，對除訓練集以外的情況進行預測，將預測得到的修正系數帶回到RANS 求解器中，提高RANS 湍流模型的預測精度.

2. 2流場反演

針對修正后的SA 湍流模型，設置一個目標函數使SA 湍流模型求解出的氣動參數與實驗測量的數據更加吻合，本文選擇升力作為優化目標，則該目標函數的表達式如式（3）.

（3）解碼器階段. 由MLP 組成的解碼器將處理階段最后一步的節點特征vi 轉換為一個或多個輸出特征，在本文中只有一個輸出，即SA 模型的修正項β.

3實驗結果與分析

3. 1針對S809翼型的網格生成及流場反演

本文選取了S809翼型作為研究對象，S809翼型是美國國家可再生能源實驗室（NREL）為水平軸風力機設計的翼型，實驗數據充分. 為了驗證所提出方法的有效性，首先生成了網格單元數為3000的一套網格，隨后遵循同一網格拓撲，分別在x 方向和y方向等比加密，構造了網格單元分別為12 000、21000、30000等3套網格. 圖5 展示了網格單元數分別為3000及30000的網格.

來流雷諾數2×106，分別在8. 2°、12. 2°、14. 2°等3 個攻角下進行流場反演. 圖6 展示了在3個攻角下計算得到的流場圖.

從圖6 可以看到，在8. 2°攻角下，由于攻角較低，并沒有發生流動分離現象，但是隨著攻角增大到12. 2°，在翼型尾部發生了流動分離，在翼型尾部產生了一個回流區. 攻角進一步增大到14. 2°，回流區進一步擴大，流動分離更加明顯.

通過在8. 2°、12. 2°、14. 2°等3 個攻角下進行流場反演，得到SA 湍流模型修正項β 的空間分布.如圖7 所示.

從圖7 可以看出，在8. 2°攻角，由于攻角較小，沒有發生流動分離的情況下，修正項β 也只在翼型近壁面附近發生變化，且其值變化的范圍也比較小. 隨著攻角增大到12. 2°，開始發生流動分離，β的變化區域隨之改變. 進一步增大攻角，流動分離現象更劇烈的情況下β 值變化的范圍也進一步增大. 總的來說，β 只是在翼型的上表面以及尾部發生了變化，這與流場分離區域大致符合；而在翼型下表面及遠場區域，β 基本沒有發生改變. 表1 總結了30 000 網格下流場反演在各個攻角下計算出的升力系數.

可以看到SA 模型難以準確預測復雜流動情況下的升力系數，通過流場反演，修正了原始的SA 湍流模型，使其能夠較為準確預測升力系數.

3. 2模型訓練及預測結果對比

通過流場反演得到了SA 模型修正項β 的空間分布，接下來需要運用機器學習算法訓練模型，模型的輸入是從流場中提取的具有物理意義的特征η，模型輸出是修正項β，訓練好的模型能夠預測其他流動情況下的修正，代入CFD 求解器中能夠改進SA 模型對于升力系數的預測. 利用在8. 2°、12. 2°、14. 2°攻角下進行流場反演得到的數據訓練模型，對9. 0°、10. 2°、12. 5°、13. 2°攻角下進行預測.在之前的研究中，大多采用MLP 預測修正項，本文使用PyTorch 深度學習框架搭建了一個8 層（包含輸入輸出層）的MLP 網絡，每層包含64 個神經元.搭建的基于Encoder-Process-Decoder 結構的圖神經網絡模型超參數下，Encoder 模塊與Decoder 模塊中都使用包含2 層32 維的MLPs 結構，Process階段由3 個GN-block 構成，每個GN-block 中的Edge-block 和Node-block 都由2 層32 維的MLP 構成. 使用相同的初始學習率0. 001，相同的優化器Adam 更新權重，相同的激活函數ReLU，訓練相同數量的5000 個epoch，采用相同的損失函數MSE.圖8 展示了在網格單元數為3000 的網格上14°攻角下預測的修正項β 的相對誤差云圖. 從圖8 展示的誤差分布可以看出，誤差主要還是集中在翼型的上表面，雖然看上去GNN 模型所預測的誤差范圍更廣，但是實際上GNN 預測的誤差卻比MLP 低一個數量級，說明融入空間信息的GNN 模型能更好地預測β，能取得更好的效果.

表2 總結了3000網格下MLP與GNN對各個攻角計算出的升力系數及其相對誤差.

從表2 可以看到，相較于MLP，GNN模型所預測的升力系數更接近實驗數據，預測的相對誤差也更小，更能夠改進SA 湍流模型對升力系數的預測.

圖9展示了網格數分別為3000、12000、21000、30000 這4套網格上的結果. 圖9 中橫坐標aoa 表示攻角大小，括號中表示在不同網格上計算結果，縱坐標表示求解出的升力系數. 可以看出在保持MLP 與GNN 相同參數量情況下，總的來說在4 套網格上GNN 求解出的升力系數都要比MLP 更加接近實驗測量出的升力系數，這是由于GNN 在預測修正項的過程中引入了空間信息，利用消息傳遞機制，收集周圍多個點的特征信息一起用來預測修正項.

3. 3針對GNN 的消息傳遞機制的改進

通過進一步研究發現輸入特征δ 的分布與β很相似，都是集中在翼型表面以及尾部附近，如圖10. 為了使GNN 能夠更加敏銳地捕捉到翼型附近β 的變化，因此構造了新的消息傳遞機制. 在新的消息傳遞機制中利用特征δ 增添了兩個權重wv_i和we_ij，新的消息傳遞更新過程如式（18）和式（19）所示.

圖9 展示了網格數分別為3000、12000、21000、30 000 這4 套網格上的結果. 圖9中橫坐標aoa表示攻角大小，括號中表示在不同網格上計算結果，縱坐標表示求解出的升力系數. 可以看出在保持MLP與GNN 相同參數量情況下，總的來說在4 套網格上GNN 求解出的升力系數都要比MLP 更加接近實驗測量出的升力系數，這是由于GNN 在預測修正項的過程中引入了空間信息，利用消息傳遞機制，收集周圍多個點的特征信息一起用來預測修正項.

在更新點上特征與邊上特征時分別乘上權重w_{v_i}、w_{e_ij}，權重的值通過特征δ 來計算，計算表達式如式（20）和式（21）

從特征δ的分布圖10可以看出在遠場δ 的取值接近于0，經過一個正切函數之后賦予一個接近于零的權重，但是在翼型表面附近，經過正切函數之后就得到了一個非零的權重，引入權重之后弱化了GNN模型在訓練過程中遠場部分的消息傳遞，突出了翼型表面附近的消息傳遞，使得GNN對壁面附近的空間信息更加敏感，進而改進了GNN模型的效果. 改進消息傳遞機制后的效果如圖11.

從圖11 可以看到，在改進GNN 的消息傳遞機制后，GNN 的效果得到進一步的改進，并且通過加權的方式改變消息傳遞機制，既沒有增大計算開銷，也沒有增大時間復雜度，是一種較理想的方式.

4結論

本文基于FIML框架，針對目前該框架內普遍采用的多層感知機網絡對于湍流空間相關性表征不足的缺陷，提出使用圖神經網絡進行湍流建模，搭建了基于Encoder-Process-Decoder結構的圖神經網絡模型. 并將網格點之間的相對位移向量以及兩點間的歐式距離等空間信息編碼成邊上的特征，通過消息傳遞對點上特征和邊上特征進行聚合. 實驗結果表明，在保證參數量與傳統MLP 網絡結構相當的同時，提升了預測的準確率. 進一步地，為了使GNN能夠更加準確地捕捉到翼型附近的修正項變化，本文改進了GNN 訓練過程的消息傳遞機制. 通過訓練策略的調整，使得在參數量相當的情況下，GNN 的效果得到進一步提升. 這為圖神經網絡的發展提供了一個新的思路：即考慮輸入特征與輸出特征的分布情況來合理設計這樣一種加權函數，突出圖中重點部分的消息傳遞，進而改進圖神經網絡的預測效果.

盡管引入空間信息的FIML 方法在S809翼型中取得了較之前更好的效果，但是本文提出方法的泛化能力仍有進一步提升的空間. 首先，訓練的模型推廣應用到不同雷諾數時預測結果與實驗結果差異較大；其次，訓練的模型應用到不同幾何外形時結果也不夠理想. 這主要是因為在不同雷諾數、不同的幾何外形情況下，邊界層內的流動特征差異較大，在給定數據集上訓練得到的修正模型，在推廣應用于其它情形時，仍有局限性，這也是后續將要重點研究的方向. 除此之外，本文的湍流模型修正系數的優化目標暫時僅考慮翼型的升力系數，實際上翼型的氣動性能還包括阻力系數，力矩系數等，而其對應的流場特性還包括壓力、速度等的分布. 要準確預測復雜湍流流動，這些物理量也是應該關注的. 因此，在后續的研究中，我們將針對更寬廣的優化目標開展傳統湍流模型的修正方法研究.