趙鄭豪, 王之歡, 朱洪強(qiáng)

(南京郵電大學(xué)理學(xué)院, 南京 210023)

非線性雙曲守恒律方程的數(shù)值解通常會(huì)出現(xiàn)間斷,導(dǎo)致數(shù)值計(jì)算中產(chǎn)生振蕩, 造成不穩(wěn)定性, 給數(shù)值模擬帶來很大的困難, 而RKDG方法(Runge-Kutta discontinuous Galerkin methods)[1]是求解這類方程的主流方法之一, 其主要過程為: 1) 通過設(shè)計(jì)壞單元指示子快速定位間斷位置, 并標(biāo)記出附近的壞單元; 2) 重構(gòu)壞單元上的數(shù)值解控制數(shù)值振蕩．RKDG方法的另一個(gè)特點(diǎn)是近似解不要求單元間的連續(xù)性, 故該方法能在各類網(wǎng)格上使用, 包括帶有懸點(diǎn)的非協(xié)調(diào)網(wǎng)格．目前已有許多關(guān)于壞單元指示子的研究工作, 在較早的工作中, Qiu等[2]數(shù)值比較了7個(gè)常用的壞單元指示子, 但并沒有得到一個(gè)對(duì)多數(shù)問題表現(xiàn)一致最優(yōu)的壞單元指示子; Dumbser等[3]對(duì)2014年之前的壞單元指示子工作做了比較詳盡的回顧．近幾年, Vuik等[4]提出一種根據(jù)不同壞單元指示變量能自動(dòng)選擇參數(shù)的壞單元指示子; Gao等[5]在此基礎(chǔ)上進(jìn)行了拓展優(yōu)化; Fu等[6]提出一個(gè)既簡(jiǎn)單又不含敏感參數(shù)的壞單元指示子; 筆者[7]則將此工作推廣到了h自適應(yīng)網(wǎng)格上．近些年, 機(jī)器學(xué)習(xí)在這一領(lǐng)域也得到了廣泛應(yīng)用, 如基于神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)的壞單元指示子[8-11]等．

間斷區(qū)域和連續(xù)區(qū)域內(nèi)單元的指示變量值往往相差較大, 筆者[12]基于這個(gè)觀察設(shè)計(jì)了一個(gè)基于K均值聚類的壞單元指示子, 該工作通過將所有單元?jiǎng)澐譃槿舾蓚€(gè)局部小模板, 讓每個(gè)模板中至少包含一個(gè)好單元, 再利用統(tǒng)計(jì)方法中常用的K均值聚類方法將模板中的單元最終分成兩類, 然后使用一個(gè)新設(shè)計(jì)的模板分類準(zhǔn)則來判斷模板中是否存在壞單元, 從而確定每一個(gè)單元是否為壞單元．這個(gè)K均值聚類壞單元指示子對(duì)局部收集來的指示變量值運(yùn)用數(shù)據(jù)分析方法做壞單元判斷, 避開了對(duì)網(wǎng)格剖分的依賴, 可以用于各種類型的網(wǎng)格, 且在數(shù)值試驗(yàn)中能準(zhǔn)確判斷間斷位置．然而,此方法含有2個(gè)實(shí)數(shù)參數(shù), 雖然它們不依賴測(cè)試問題, 但參數(shù)的確定需要大量的數(shù)值測(cè)試, 影響了該方法的實(shí)際可用性．為了解決這個(gè)問題, 本文擬基于K均值聚類壞單元指示子[12]的框架, 在聚類之前首先對(duì)指示子變量值進(jìn)行歸一化來消除不同算例中指示變量取值范圍不同的影響,然后利用極差判斷模板中的單元是否存在壞單元．這種做法不僅簡(jiǎn)單、高效,而且僅含1個(gè)參數(shù), 能極大地降低確定參數(shù)的復(fù)雜度, 有利于該方法在實(shí)際工作中的應(yīng)用．

1 RKDG方法

2 改進(jìn)的K均值聚類壞單元指示子

受文獻(xiàn)[12]壞單元指示子框架的啟發(fā), 本文將壞單元指示子仍然分為兩部分: 1) 為了讓探測(cè)壞單元在解結(jié)構(gòu)更為簡(jiǎn)單的局部區(qū)域進(jìn)行, 將計(jì)算單元按一定規(guī)則劃分為一些模板; 2) 對(duì)每個(gè)模板進(jìn)行壞單元探測(cè)．具體步驟:

i) 對(duì)于第(1)部分, 模板劃分的方法參見文獻(xiàn)[12]．因?yàn)槿呛脝卧腿菈膯卧０迳现甘咀兞康闹翟跀?shù)據(jù)表現(xiàn)上不易區(qū)分, 故要求每個(gè)模板中至少有1個(gè)好單元．同時(shí), 為了節(jié)省計(jì)算時(shí)間, 模板大小(即模板中包含的單元數(shù))越小越好, 故一維時(shí)模板尺寸取為5, 二維矩形網(wǎng)格時(shí)模板尺寸取為4×4[12]．

ii) 對(duì)于第(2)部分, 需要對(duì)每個(gè)單元判斷是否為壞單元．這部分又可細(xì)分為: 歸一化處理、模板篩選和壞單元篩選．由于壞單元指示變量對(duì)于不同的算例或不同的DG多項(xiàng)式次數(shù)k, 其取值范圍各不同, 導(dǎo)致壞單元指示子中的參數(shù)可能會(huì)依賴于算例或k值, 給數(shù)值應(yīng)用帶來麻煩．受文獻(xiàn)[13]的啟發(fā), 首先對(duì)指示變量的值進(jìn)行歸一化處理．假設(shè)計(jì)算單元總數(shù)為n, 給定指示變量I, 所有單元上的指示變量值分別記作I1,I2,…,In．定義Imax=max{I1,I2,…,In},Imin=min{I1,I2,…,In}．對(duì)每個(gè)單元的指示變量值Ii做線性變換

Ti=(Ii-Imin)/(Imax-Imin)

,

(1)

指示變量須進(jìn)行歸一化處理, 否則對(duì)不同情況下參數(shù)σ的最優(yōu)取值明顯不同[13], 這給實(shí)際應(yīng)用帶來不便．經(jīng)歸一化處理的壞單元指示子可使用多種壞單元指示變量, 本文采用單元邊界跳量指示變量[14]．假設(shè)目標(biāo)單元為K0, 且有G個(gè)直接相鄰單元, 分別是K1,K2,…,KG．若pj(j=0,1,…,G)是單元Kj上的DG近似解, 則單元邊界跳量指示變量

注與文獻(xiàn)[12]的做法相比, 本文算法中去掉了原來含有2個(gè)參數(shù)的過濾步驟和模板分類判斷步驟, 改為使用歸一化來優(yōu)化數(shù)據(jù), 結(jié)合只含有1個(gè)參數(shù)的極差法做模板分類判斷, 能簡(jiǎn)單、高效地篩選出含有壞單元的模板．由于新的壞單元指示子只涉及1個(gè)參數(shù), 少量的測(cè)試就能確定參數(shù)取值, 極大地方便了該法的數(shù)值應(yīng)用．

3 數(shù)值結(jié)果

以歐拉方程組作為模型方程, 采用經(jīng)典的算例來展示改進(jìn)的壞單元指示子的效果．若方程在二維情形時(shí)的形式為

其中ρ為密度,E為總能量,p為壓力,u和v分別為x和y方向上的速度分量, 這些量滿足p=(γ-1)[E-0.5ρ(u2+v2)], 式中γ=1.4．改進(jìn)的壞單元指示子基于密度和能量分別探測(cè)壞單元, 一個(gè)單元只要有一次被探測(cè)為壞單元, 最終就會(huì)被標(biāo)記為壞單元．

例1Lax問題．有初值條件

取網(wǎng)格單元大小為Δx=1/40, 圖1(a)和(b)為t=1.3時(shí)的密度數(shù)值解, 圖1(c)和(d)為不同時(shí)刻壞單元的位置圖(壞單元?dú)v史圖)．計(jì)算結(jié)果表明, 改進(jìn)的壞單元檢測(cè)方法能準(zhǔn)確標(biāo)記出間斷的位置, 數(shù)值解沒有產(chǎn)生振蕩．

圖1 Lax問題在t=1.3時(shí)的結(jié)果Fig.1 Results at t=1.3 for the Lax problem

例2爆炸波問題．初值條件為

圖2給出了Δx=1/800的網(wǎng)格在t=0.038時(shí)的數(shù)值解．從圖2可以看到, 改進(jìn)的壞單元指示子準(zhǔn)確探測(cè)到了2個(gè)大間斷, 且數(shù)值解無振蕩．

圖2 爆炸波問題在t=0.038時(shí)的結(jié)果Fig.2 Results at t=0.038 for the blast-wave problem

例3雙馬赫反射問題．設(shè)計(jì)算區(qū)域?yàn)閇0,4]×[0,1], 底邊1/6≤x≤4區(qū)域處為固壁．當(dāng)過點(diǎn)(1/6,0)且與固壁成60°的激波以10馬赫的速度沿垂直方向撞擊固壁時(shí), 設(shè)激波前方?jīng)]有受到擾動(dòng)的空氣密度為1.4, 壓力為1; 底側(cè)的固壁區(qū)域使用反射邊界條件, 其他區(qū)域使用精確的波后邊界條件; 頂部邊界使用10馬赫激波運(yùn)動(dòng)的精確邊界條件, 左右邊界依次使用入流邊界條件和出流邊界條件．圖3給出了t=0.2,Δx=Δy=1/240時(shí)得到的密度等高線和壞單元．結(jié)果顯示, 雖然在k=1的情形中間斷位置沒有被完全探測(cè)出來, 但所有的數(shù)值解沒有振蕩．改進(jìn)的壞單元指示子在確保數(shù)值解沒有振蕩的前提下,標(biāo)記了盡可能少的壞單元．

圖3 雙馬赫反射問題在t=0.2時(shí)的結(jié)果Fig.3 Results of the double Mach refection problem at t=0.2

例4前臺(tái)階問題．問題起始于一個(gè)長(zhǎng)度為3, 高度為1的風(fēng)洞中, 風(fēng)洞底部有一高為0.2、離風(fēng)洞左邊端口0.6并一直延伸到風(fēng)洞右端口的臺(tái)階; 風(fēng)洞壁處使用反射邊界條件, 左右邊界分別使用入流和出流邊界條件．圖4給出了t=4.0, Δx=Δy=1/160時(shí)得到的密度等高線以及壞單元．從圖4可以看到, 壞單元探測(cè)準(zhǔn)確, 數(shù)值解并沒有振蕩．

圖4 前臺(tái)階問題在t=0.4時(shí)的結(jié)果Fig.4 Results of the forward-facing step problem at t=0.4

為了驗(yàn)證改進(jìn)的方法仍然適用于不同的壞單元指示變量, 圖5給出這個(gè)算例使用Fu-Shu指示變量[14]的計(jì)算結(jié)果．從圖5可以看到, 壞單元檢測(cè)準(zhǔn)確, 數(shù)值解沒有振蕩．

為了說明改進(jìn)的壞單元指示子在不同粗細(xì)網(wǎng)格上的表現(xiàn), 表1列出了這個(gè)算例2個(gè)指示變量在兩套網(wǎng)格上的壞單元平均百分比．結(jié)果顯示, 當(dāng)網(wǎng)格尺寸減半時(shí), 壞單元平均百分比在所有情況下都大約減小到原來的一半, 說明改進(jìn)的壞單元指示子在不同粗細(xì)網(wǎng)格上的表現(xiàn)良好．

表1 前臺(tái)階問題壞單元平均百分比