稀疏發(fā)現(xiàn)湍流邊界層壁面函數(shù)公式

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

(1. 西南科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,四川綿陽621000;2.中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算空氣動(dòng)力研究所,四川綿陽621000)

1 引言

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

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

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

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

2 壁面函數(shù)

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

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

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

(1)

其中

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

但計(jì)算流體學(xué)界對(duì)統(tǒng)一壁面函數(shù)這一經(jīng)驗(yàn)公式的正確性一直存在爭(zhēng)議,本文從實(shí)驗(yàn)數(shù)據(jù)中獲取方程以驗(yàn)證該方程的準(zhǔn)確性。除了溫度壁面函數(shù),工程上還會(huì)用到速度壁面函數(shù),本文主要研究溫度壁面函數(shù),有關(guān)速度壁面函數(shù)以及標(biāo)準(zhǔn)壁面函數(shù)等其它內(nèi)容可以參考文獻(xiàn)21到23。

3 稀疏辨識(shí)算法

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

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

y=f(x)

(2)

其中f代表一個(gè)非線性函數(shù),根據(jù)收集到的數(shù)據(jù)可以得到以下矩陣關(guān)系

(3)

構(gòu)造一個(gè)n列的候選非線性函數(shù)的函數(shù)庫Μ(x)如式(4)所示。Μ(x)可以由多種不同的函數(shù)項(xiàng)構(gòu)成,候選庫的設(shè)計(jì)根據(jù)領(lǐng)域知識(shí)由領(lǐng)域?qū)＜掖_定。

(4)

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

y=Μ(x)η

(5)

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

(6)

若η中某一特定項(xiàng)的系數(shù)不為零,則其相對(duì)應(yīng)的函數(shù)項(xiàng)在本模型中是活躍的。函數(shù)f通常只包含幾個(gè)重要的項(xiàng),使得它在可能的函數(shù)空間中表示出強(qiáng)稀疏性。

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

(7)

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

(8)

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

4 壁面數(shù)據(jù)挖掘

4.1 實(shí)驗(yàn)準(zhǔn)備

本文所用數(shù)據(jù)來自于壁面網(wǎng)格實(shí)驗(yàn),整個(gè)網(wǎng)格實(shí)驗(yàn)在Ubuntu16.04系統(tǒng)下結(jié)合OpenFOAM5.0版本軟件進(jìn)行,對(duì)5種不同馬赫數(shù)的平板邊界層進(jìn)行了數(shù)值模擬,詳細(xì)的流動(dòng)情況見表1。本工作針對(duì)高馬赫數(shù)下壁面函數(shù)的挖掘,因此選取的均為馬赫數(shù)大于5的算例,且表中算例是平板邊界層中的經(jīng)典算例[24]。馬赫數(shù)的增加會(huì)增強(qiáng)湍流結(jié)構(gòu)的變化,不同的馬赫數(shù)算例可以得到多樣的數(shù)據(jù),這使得挖掘出的壁面函數(shù)更具有普適性。

表1 算例信息

Crocco-Busemann溫度統(tǒng)一壁面函數(shù)與實(shí)驗(yàn)數(shù)據(jù)存在較大差異,其結(jié)果如圖1所示。由圖1可以看出,Crocco-Busemann函數(shù)并不適合用于工程應(yīng)用獲得溫度信息。為了獲得適用于工程的溫度壁面函數(shù),對(duì)溫度實(shí)驗(yàn)數(shù)據(jù)進(jìn)行數(shù)據(jù)挖掘。

圖1 Crocco-Busemann函數(shù)與實(shí)驗(yàn)值對(duì)比圖

4.2 實(shí)驗(yàn)數(shù)據(jù)挖掘

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

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

(9)

結(jié)合式(9)的函數(shù)形式對(duì)壁面實(shí)驗(yàn)數(shù)據(jù)進(jìn)行挖掘得到系數(shù)A和B,進(jìn)而得到適用于工程的統(tǒng)一壁面函數(shù)。

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

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

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

(10)

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

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

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

(11)

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

圖2 控制A得到的統(tǒng)一函數(shù)效果對(duì)比

控制變量B取原函數(shù)值,而對(duì)修正系數(shù)A進(jìn)行計(jì)算來獲取統(tǒng)一壁面函數(shù),所得到地函數(shù)形式為

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

(12)

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

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

使用STEN和LASSO算法對(duì)表3的結(jié)果進(jìn)行回歸分析,得到修正系數(shù)A的統(tǒng)一形式:

(13)

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

圖3 控制B得到的統(tǒng)一函數(shù)效果對(duì)比

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

4.3 誤差分析及結(jié)果驗(yàn)證

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

表4 均方誤差結(jié)果表

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

(14)

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

(15)

將該結(jié)果在各算例上進(jìn)行驗(yàn)證,對(duì)各個(gè)算例的效果進(jìn)行展示,如圖4所示。

圖4 所得函數(shù)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比效果圖

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

5 總結(jié)

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

在接下來的工作中將進(jìn)一步分析挖掘出的統(tǒng)一壁面函數(shù)背后的物理機(jī)制,為該公式的工程應(yīng)用提供更多的理論支撐。