張 暉,葉 濤,王新光,萬 釗

(1. 西南科技大學計算機科學與技術學院,四川綿陽621000;2.中國空氣動力研究與發(fā)展中心計算空氣動力研究所,四川綿陽621000)

1 引言

計算流體力學技術已經在飛行器設計、流體機械設計等領域廣泛應用,能夠實現復雜外形高超聲速流動的數值模擬,成為風洞實驗的重要補充。數值模擬技術若要準確地獲得壁面處的物理信息分布,需要在壁面邊界層中布置非常細密的網格,這樣會極大地增加迭代收斂的計算步數,使迭代計算的穩(wěn)定性下降,同時還會帶來嚴重的數值剛性問題。工程上一般采用壁面函數的方法來獲得壁面熱流和速度等物理信息,壁面函數的使用在很大程度上放寬了近壁面網格密度,降低了對網格的依賴性,加速了收斂過程,在處理邊界層計算問題時,壁面函數發(fā)揮著重要作用[1]。

邊界層沿壁面垂直方向可以分成三個子層:粘性底層,過渡層,對數律層,不同層具有不同的流動狀態(tài)。由于過渡層的復雜性,標準壁面函數通常設計為兩層模型,分別給出了粘性底層和對數律層的壁面規(guī)律,并未對過渡層做出任何處理,常常會產生速度和溫度的不連續(xù)。而統一壁面函數與實驗值相差較大,無法準確地得到壁面附近的物理量變化情況,僅適用于理論研究中觀察其趨勢,并不適用于工程應用。

機器學習能夠發(fā)現數據中存在的規(guī)律,近年來,數據驅動發(fā)現控制方程的機器學習方法引起了廣泛的關注。由于計算量較小,在多種機器學習算法中,稀疏回歸方法在發(fā)現各種問題的控制方程方面表現出巨大的潛力。使用稀疏回歸可以從預先定義的大候選庫中識別構成控制方程的少量項,得到一個簡約模型。Brunton[2]提出了一個非線性動力學稀疏辨識框架來發(fā)現動態(tài)系統的控制方程。Rudy[3]在Brunton的基礎上提出了順序閾值嶺回歸算法,產生了一個非線性動力學的偏微分方程函數辨識框架。在這之后,大量的文獻開始研究稀疏回歸在控制方程數據驅動發(fā)現中的應用[4-13]。本工作之前的成果表明稀疏回歸能夠用于仿真數據的控制方程獲取并從中得到壁面函數[14]。除了稀疏回歸方法,其它技術也被用于執(zhí)行數據驅動的控制方程的發(fā)現,如高斯過程和神經網絡[15-17]。Maziar[18]等提出了一個利用高斯過程來發(fā)現控制方程的框架。Long[19]等在其之前的工作基礎[20]上提出了一種數字和符號混合的深度神經網絡,通過學習復雜系統的動力學來揭示隱藏的偏微分方程模型。這些方法方法不會對偏微分方程中的無關項進行懲罰,因此不能產生簡約的模型。

本文基于彈性網絡法提出了一種新的稀疏回歸算法,稱為順序閾值彈性網絡(Sequentially Threshold Elastic Net, STEN),并首次將稀疏回歸方法用于壁面實驗數據識別并得到了適用于工程的壁面函數。STEN結合了LASSO和嶺回歸的優(yōu)勢,在保證穩(wěn)定性的同時也具有可快速計算的優(yōu)點。將提出的STEN算法和傳統的稀疏回歸LASSO算法分別用于壁面實驗數據擬合。計算結果表明,LASSO和STEN兩種算法在壁面數據擬合中均表現出色,且STEN算法比LASSO算法相比結果更準確,效果更好。

2 壁面函數

湍流壁面函數是一種使用預定函數作為邊界條件來解決近壁面復雜現象的方法。壁面函數允許在壁面附近使用相對粗糙的網格單元,從而節(jié)省計算時間。應用簡化的雷諾平均Navier-Stokes方程,可以將不同變量(速度、溫度等)表示為關于壁面距離的函數,并將邊界層中的流動信息作為第一個節(jié)點的邊界條件來進行數值模擬,這些函數稱為壁面函數。

邊界層不同層的分子粘性和湍流粘性作用不同,其中過渡層分子粘性和湍流粘性的作用相當,流動狀態(tài)比較復雜,很難用一個公式或定律描述,且其厚度較小,工程上一般忽略過渡層的影響。所以,標準壁面函數常設計為兩層模型。采用兩層模型,易造成物理信息分布的不連續(xù)性,因此有學者提出了統一壁面函數。Crocco-Busemann提出了在對數層和粘性底層內統一有效的可壓縮流動邊界層壁面規(guī)律,其中溫度統一壁面函數如式(1)所示[21-23]

T/Tw=1+βu+-Γ(u+)2

(1)

其中

u+為無量綱速度,μw為壁面粘性系數,qw為壁面熱流,ρw為壁面密度,kw為壁面熱傳導系數,uτ為摩擦速度,Cp為壓力系數,r為壁面律常數,Tw為壁面溫度,T為溫度。

但計算流體學界對統一壁面函數這一經驗公式的正確性一直存在爭議,本文從實驗數據中獲取方程以驗證該方程的準確性。除了溫度壁面函數,工程上還會用到速度壁面函數,本文主要研究溫度壁面函數,有關速度壁面函數以及標準壁面函數等其它內容可以參考文獻21到23。

3 稀疏辨識算法

本文采用稀疏辨識算法來處理回歸問題。根據收集到的數據設計回歸模型的函數候選庫,并使用稀疏辨識方法從數據中發(fā)現方程。稀疏辨識很好地避免了過擬合,同時還能通過壓縮系數來降低函數的復雜度。

假設自變量和因變量間存在如下規(guī)律

y=f(x)

(2)

其中f代表一個非線性函數,根據收集到的數據可以得到以下矩陣關系

(3)

構造一個n列的候選非線性函數的函數庫Μ(x)如式(4)所示。Μ(x)可以由多種不同的函數項構成,候選庫的設計根據領域知識由領域專家確定。

(4)

Μ(x)每一列代表式(2)右邊的一個候選函數項,在許多物理系統中,這些非線性函數項只有少數在f中是活躍的,因此可以通過稀疏回歸來確定系數η=[η1,η2,η3,…,ηm]的稀疏矢量

y=Μ(x)η

(5)

η中的每一列代表一個稀疏向量。其中η=[η1,η2,η3,…,ηm]是大小為m的系數向量,m是庫Μ(x)中的候選項個數。將根據自變量數據x設計的候選函數庫Μ(x)和收集到的y帶入到(5)中,通過求解系數矩陣η來獲得函數形式

(6)

若η中某一特定項的系數不為零,則其相對應的函數項在本模型中是活躍的。函數f通常只包含幾個重要的項,使得它在可能的函數空間中表示出強稀疏性。

針對稀疏求解獲得η,可采用最小二乘法,嶺回歸、LASSO算法和順序閾值最小二乘法。最小二乘法會辨識所有函數項,使得所有系數均為非零值,可能會導致過擬合;LASSO算法壓縮系數能力較強,但穩(wěn)定性有所減弱;嶺回歸穩(wěn)定性高但壓縮系數的能力較低。而彈性網絡算法結合了LASSO和嶺回歸兩種算法的優(yōu)點,既能夠壓縮系數也能夠保證穩(wěn)定性。本文基于彈性網絡算法提出了一種新的稀疏求解算法,稱為順序閾值彈性網絡法(Sequentially Threshold Elastic Net, STEN)來求解稀疏向量。該方法是以縮小變量集為思想的壓縮估計方法,它通過設計一個懲罰函數來壓縮變量的系數并使部分函數項回歸系數為零,以此來進行變量選擇。STEN的代價函數為

(7)

其中λ1和λ2分別對L1和L2懲罰函數的權重進行正則化。方程(7)中的λ1和λ2通過對懲罰項施加更多的權重來控制稀疏性的數量,從而導致系數的收縮和避免過擬合?？梢允褂靡粋€正則化權重參數α和混合參數λ來修正方程(7)

(8)

利用λ參數控制L1和L2的凸組合,當λ=1時為LASSO算法,而λ=0為嶺回歸算法。正則化的權重α與懲罰項相乘,可以決定模型所需的稀疏強度,α=0代表一般的最小二乘法,α越大則稀疏性越強。為保證結果的稀疏性,對方程(8)結果添加一個閾值,對小于閾值的系數設為0,再將剩余系數進行遞歸計算直到非零系數收斂。該算法計算效率高,能夠通過少量迭代快速收斂得到稀疏解。

4 壁面數據挖掘

4.1 實驗準備

本文所用數據來自于壁面網格實驗,整個網格實驗在Ubuntu16.04系統下結合OpenFOAM5.0版本軟件進行,對5種不同馬赫數的平板邊界層進行了數值模擬,詳細的流動情況見表1。本工作針對高馬赫數下壁面函數的挖掘,因此選取的均為馬赫數大于5的算例,且表中算例是平板邊界層中的經典算例[24]。馬赫數的增加會增強湍流結構的變化,不同的馬赫數算例可以得到多樣的數據,這使得挖掘出的壁面函數更具有普適性。

表1 算例信息

Crocco-Busemann溫度統一壁面函數與實驗數據存在較大差異,其結果如圖1所示。由圖1可以看出,Crocco-Busemann函數并不適合用于工程應用獲得溫度信息。為了獲得適用于工程的溫度壁面函數,對溫度實驗數據進行數據挖掘。

圖1 Crocco-Busemann函數與實驗值對比圖

4.2 實驗數據挖掘

Crocco-Busemann溫度壁面函數計算結果雖然與實驗數據存在較大差異,但其整體趨勢是符合實驗結果的,同時相關理論研究也肯定了其價值。本文參考Crocco-Busemann函數形式設計了新的函數形式,如式(9)所示。在無量綱速度為零的時候,流體還未進行流動,此時壁面邊界層的溫度就是壁面本身的溫度,壁面邊界層不具有熱量梯度不會產生熱傳遞。因此在無量綱速度為零時,邊界層的溫度與壁面溫度一致,基于此理由,將常數項設置為1,與理論公式一致。

T/Tw=1+B*βu++A*Γ(u+)2

(9)

結合式(9)的函數形式對壁面實驗數據進行挖掘得到系數A和B,進而得到適用于工程的統一壁面函數。

實驗過程中結合序列最小優(yōu)化思想(Sequential minimal optimization, SMO)對STEN和LASSO算法進行優(yōu)化,將初始問題不斷地分解為子問題,再對子問題進行求解分析,最終使得全部變量都滿足條件。根據SMO算法思想對方程(9)進行控制變量處理,控制A和B其中一個值為定值,對另外一個值進行回歸計算。

首先控制修正系數A為定值,定值的選擇上延續(xù)原函數的取值,原函數中A的值為-1,方程表示如下

T/Tw=1+B*βu+-Γ(u+)2

(10)

使用平方誤差最小法計算得到各算例的B值如表2所示。

表2 控制修正系數A后各算例的B值

在流體力學中,各個物理量之間都存在可以用函數來表示的關系,而各個物理量都可以通過基本物理量來計算得到。因此選取壁面函數研究中的基本物理量作為自變量,分別是:馬赫數Ma,壁面剪切應力τw,壁面密度ρw和壁面溫度Tw。根據自變量設計了88項函數項來組成候選函數庫,庫中包含常數項、一次項、二次項、三次項、平方根項,三角函數項等等。函數侯選庫設計完成后,使用STEN和LASSO算法進行回歸分析,STEN算法在正則化權重參數α=0.09和混合參數λ=0.8時可以得到統一結果

(11)

把回歸得到的B的統一形式帶入到方程(10)中進行計算,并將計算結果與實驗結果和理論函數進行比較,部分算例的對比圖如圖2所示,圖中只展示了Case3算例的效果圖,其余算例效果與此相似。由圖2可知,兩種稀疏算法在控制A情況下所得的統一函數效果均不好。其誤差雖然較原理論函數有所降低,但與實驗結果之間的誤差仍然較大,并不能很好地匹配實驗值。

圖2 控制A得到的統一函數效果對比

控制變量B取原函數值,而對修正系數A進行計算來獲取統一壁面函數,所得到地函數形式為

T/Tw=1+βu++A*Γ(u+)2

(12)

原函數中B值為常數1,選取與原函數一致的B值,能夠得到如式(12)表示的函數表達式。使用平方誤差最小法計算得到各算例的A值如表3所示

表3 控制修正系數B后各算例的A值

使用STEN和LASSO算法對表3的結果進行回歸分析,得到修正系數A的統一形式:

(13)

將所得A值帶入到(12)中進行計算并與實驗數據進行對比。由圖3可知,兩種算法的結果均與實驗數據匹配度較高,相對于原溫度壁面函數效果均有所提升。圖3僅展示了Case3算例的效果圖,其余算例效果與此相似。

圖3 控制B得到的統一函數效果對比

從所得結果可知,B取理論公式值,計算控制變量A所得結果對原壁面函數改進最為明顯,與實驗數據的匹配度較高。進一步對控制變量A所得結果進行誤差分析,選取誤差最小的作為最終回歸結果并進行驗證。

4.3 誤差分析及結果驗證

本文的工作是從壁面數據中抽取出一個統一的壁面函數,更關注壁面函數的計算結果與真實值數據之間的整體匹配度。因此對壁面實驗數據使用不同稀疏算法得到的統一壁面函數進行誤差分析,使用均方誤差來衡量誤差大小。誤差計算結果如表4所示,雖然STEN和LASSO算法所得結果與實驗值誤差均較小,但相比之下,STEN算法誤差更小,效果較LASSO算法更好。

表4 均方誤差結果表

根據誤差分析結果選取誤差更小的STEN算法所得結果為最終結果。雖然式(13)的A表達式有表達出一定的物理意義,若0.0013這樣的小數能用物理閉合系數來表示則能夠更好的體現出其與物理量的相關性。其中常數-1.0021進行四舍五入為-1,對0.0013系數使用壁面常見閉合系數進行表示可以得到,常見的物理閉合系數有κ,B,Cp,r,Prl,Prt等。

(14)

將A的最終表達式(14)帶入到(12)中得到最終確定的統一壁面函數形式為式(15)。

(15)

將該結果在各算例上進行驗證,對各個算例的效果進行展示,如圖4所示。

圖4 所得函數與實驗數據對比效果圖

式(15)為本文根據原壁面函數形式結合SMO算法優(yōu)化后STEN回歸算法對實驗數據,進行回歸分析得到的最終函數形式。所得結果不僅能夠較好地匹配實驗數據,使得后期工程應用中關于壁面的處理成本降低,同時也能提高數值模擬的效率。

5 總結

本文提出了一種基于彈性網絡的稀疏回歸算法STEN,并將其對壁面實驗數據進行數據挖掘得到控制方程。STEN與傳統的LASSO稀疏回歸算法相比,二者在壁面實驗數據的挖掘中效果均較好。誤差分析結果表明,STEN算法較LASSO算法誤差更小,準確度更高,STEN算法對于壁面數據更具適用性。通過對實驗數據進行擬合得到了統一壁面函數,其與實驗數據匹配度高,對比理論壁面函數效果有較大提升。將本文研究所得的統一壁面函數應用于工程應用中可以大大縮短數值模擬的時間,提高工程效率。

在接下來的工作中將進一步分析挖掘出的統一壁面函數背后的物理機制,為該公式的工程應用提供更多的理論支撐。