非穩(wěn)態(tài)Hamilton-Jacobi方程的7階加權(quán)緊致非線性格式1)

胡迎港 蔣艷群 黃曉倩

(西南科技大學(xué)數(shù)理學(xué)院,四川綿陽(yáng) 621010)

引言

d維空間下非穩(wěn)態(tài)Hamilton-Jacobi (HJ)方程定義如下

其中,t>0 ,x=(x1,x2,···xd)T,這類方程在物理學(xué)、流體力學(xué)、圖像處理、微分幾何、金融數(shù)學(xué)、最優(yōu)化控制理論等諸多領(lǐng)域中有著廣泛應(yīng)用[1-4].其解的性質(zhì)很特殊,如弱解存在但不唯一,即使初值和Hamilton 函數(shù)充分光滑,方程解的導(dǎo)數(shù)也可能出現(xiàn)間斷[5].Crandall 等[6]首次提出了HJ 方程的黏性解的概念,并證明了黏性解的存在唯一性.HJ 方程的性質(zhì)與雙曲守恒律方程十分相似,因此,可以將雙曲守恒律方程的一些成熟數(shù)值方法[7-9]推廣用于求解HJ 方程.

目前已有較多的高精度算法被推廣應(yīng)用于HJ 方程.如Osher 等[10]通過(guò)選擇最光滑的候選模板重構(gòu)數(shù)值通量方法設(shè)計(jì)了HJ 方程的二階本質(zhì)無(wú)振蕩 (ENO) 格式.Jiang 等[11]通過(guò)將所有候選模板的重構(gòu)通量進(jìn)行了非線性加權(quán)得到了HJ 方程的5 階加權(quán)本質(zhì)無(wú)振蕩 (WENO) 格式.相比于ENO 格式,WENO 格式在光滑區(qū)域能達(dá)到更高階精度,同時(shí)在間斷區(qū)域恢復(fù)為ENO 格式.Bryson等[12]設(shè)計(jì)了HJ 方程的5 階中心加權(quán)本質(zhì)無(wú)振蕩(CWENO)格式.Qiu 等[13]提出了HJ 方程的基于Hermite 多項(xiàng)式的5 階HWENO 格式,該格式在重構(gòu)通量時(shí)具有更好的緊致性.Ha 等[14]為改善格式的間斷捕捉能力,采用拉格朗日型指數(shù)多項(xiàng)式建立了一種新的6 階WENO格式.相比于其他類型的WENO 格式,該格式具有更低的計(jì)算成本.Cheng 等[15]基于經(jīng)典的5 階WENO 格式,采用4 個(gè)二次多項(xiàng)式的非線性凸組合重構(gòu)數(shù)值通量,提高了格式的精度 (在光滑區(qū)域能達(dá)到6 階) 和健壯性.Zhu 等[16]基于原始的5 階HWENO 格式,通過(guò)采用大小不同的候選模板構(gòu)造了HJ 方程的更加簡(jiǎn)單有效的HWENO 格式.Rathan[17]提出了一種新的5 階WENO 格式求解HJ 方程,該格式通過(guò)設(shè)計(jì)L1范數(shù)型的非線性權(quán)重來(lái)提高精度和分辨率.Abedian[18]采用對(duì)稱的5 階WENO 格式求解HJ 方程.Kim 等[19]基于指數(shù)多項(xiàng)式構(gòu)造了一種新的3 階WENO 格式.其特點(diǎn)是可以很好地分辨出光滑區(qū)域和間斷區(qū)域.

高階精度加權(quán)緊致非線性格式 (WCNS) 是另一種求解雙曲守恒律方程的有效方法.Deng 等[20]基于WENO 格式的構(gòu)造思想,在緊致非線性格式[21](CNS) 中引入非線性WENO 插值技術(shù),提出了高階精度WCNS 格式.該格式被證實(shí)具有高分辨率和強(qiáng)間斷捕捉能力等良好特性.Nonomura 等[22]和Zhang 等[23]分別對(duì)WCNS 格式進(jìn)行了改進(jìn),提出7 階和9 階精度的WCNS 格式.并通過(guò)數(shù)值試驗(yàn)證明了格式的高分辨率和強(qiáng)間斷捕捉能力.為了進(jìn)一步提高WCNS 格式在復(fù)雜流體計(jì)算中的應(yīng)用,Deng等[24]構(gòu)造了混合半節(jié)點(diǎn)和節(jié)點(diǎn)的加權(quán)緊致非線性格式 (HWCNS),使得格式在間斷區(qū)域有了更高的分辨率.Nonomura 等[25]設(shè)計(jì)一種更加健壯的混合型WCNS 格式.Wong 等[26]提出局部耗散加權(quán)緊致格式 (WCNS-LD).該格式在光滑區(qū)域使用耗散較小的中心插值模板,在包含間斷區(qū)域使用耗散較大的迎風(fēng)插值模板以抑制非物理振蕩.Zhao 等[27]提出一種新的可調(diào)參數(shù)的光滑度量指標(biāo),進(jìn)一步提高了WCNS 格式的間斷捕捉能力和健壯性.洪正等[28]采用5 階WCNS 格式對(duì)各向同性湍流通過(guò)正激波的情形進(jìn)行直接數(shù)值模擬.Jiang 等[29]為HJ 方程設(shè)計(jì)了5 階精度的WCNS 格式,該格式相比于同階WENO 格式具有更好的收斂性和分辨率.Hiejima[30]基于TENO 插值思想,建立WCNS-T 格式.相比傳統(tǒng)的WCNS 格式,WCNS-T 格式在強(qiáng)不連續(xù)和高頻間斷的捕捉上更具優(yōu)勢(shì).Jiang 等[31]針對(duì)等熵Navier-Stokes 方程組提出半隱式時(shí)間推進(jìn)的WCNS格式,避免顯式WCNS 格式嚴(yán)格的CFL 穩(wěn)定性限制.Ma 等[32]構(gòu)造了非線性緊致插值的WCNS 格式,數(shù)值驗(yàn)證新的WCNS 格式可用于大渦模擬.

本文在Jiang 等[29]提出的5 階精度WCNS格式基礎(chǔ)上,進(jìn)一步構(gòu)造了非穩(wěn)態(tài)HJ 方程的7 階精度WCNS 格式.HJ 方程的Hamilton 數(shù)值通量采用具有單調(diào)性的Lax-Friedrichs 方法計(jì)算,一階空間導(dǎo)數(shù)的左、右極限值采用高階精度混合節(jié)點(diǎn)和半節(jié)點(diǎn)型的8 階中心差分格式計(jì)算,半節(jié)點(diǎn)函數(shù)值采用7 階WENO 型非線性插值方法計(jì)算.3 階TVD Runge-Kutta 方法[33]用于HJ 方程的時(shí)間離散.最后通過(guò)一系列一維和二維數(shù)值算例驗(yàn)證WCNS 格式在精度,分辨率和間斷捕捉能力等方面的特性.

1 數(shù)值方法

1.1 一維HJ 方程的WCNS 格式

考慮如下形式的一維非穩(wěn)態(tài)HJ 方程

為簡(jiǎn)單起見(jiàn),采用均勻網(wǎng)格,即網(wǎng)格節(jié)點(diǎn)設(shè)置為xi=a+i?x(i=0,1,···,N).其 中,?x=(b?a)/N為 網(wǎng)格尺寸,N為網(wǎng)格數(shù).方程(2)的半離散格式為

對(duì)式(10)右端在半節(jié)點(diǎn)xi+1/2處進(jìn)行Taylor 展開(kāi),得到

其中,B k(k=0,1,2,3) 為與 ?x無(wú)關(guān)的常系數(shù).半節(jié)點(diǎn)函數(shù)值 ?i+1/2的高階線性逼近式(7)可由式(10)中4 個(gè)低階線性逼近的線性組合得到

其中,d0=1/64,d1=21/64,d2=35/64,d3=7/64 為理想線性權(quán)重.

當(dāng)所有子模板都處于光滑區(qū)域時(shí),非線性權(quán)重等于線性權(quán)重,得到7 階精度的半節(jié)點(diǎn)函數(shù)逼近;當(dāng)某些子模板包含間斷區(qū)域時(shí),令其非線性權(quán)重趨于零,防止插值穿過(guò)間斷區(qū)域,最后得到4 階精度的半節(jié)點(diǎn)函數(shù)逼近.本文選擇WENO-Z 型非線性權(quán)重,定義如下

式中,k=0,1,2,3 ,參數(shù) ε 取值為小量1 0?20,以避免分母為0.為了使得所設(shè)計(jì)WCNS 格式在光滑區(qū)域能達(dá)到7 階精度,全局光滑度量指標(biāo) τ 定義為τ =|IS3?IS2?IS1+IS0|.其中,ISk是子模板Sk(k=0,1,2,3)的光滑度量指標(biāo),定義如下

對(duì)于半離散格式(3),采用如下3 階顯式TVD Runge-Kutta 方法[33]求解

1.2 二維HJ 方程的WCNS 格式

考慮如下形式的二維非穩(wěn)態(tài)HJ 方程

這里同樣考慮均勻網(wǎng)格,即網(wǎng)格節(jié)點(diǎn)設(shè)置為xi=i?x(i=0,1,···,M),yj=j?y(j=0,1,···,N).其中,x和y方向的網(wǎng)格尺寸分別為 ?x=(b?a)/M,?y=(d?c)/N,網(wǎng)格數(shù)為M×N.

方程(20)的半離散格式為

2 WCNS 格式的精度分析

基于一維HJ 方程,式(13)改寫(xiě)為如下形式

將其代入式(5)得到

因此,WCNS 格式達(dá)到最佳7 階精度的充分條件為

3 數(shù)值實(shí)驗(yàn)

本節(jié)數(shù)值測(cè)試7 階WCNS 格式的性能.在精度測(cè)試算例中,?t=(?x)7/3,L1和L∞范數(shù)數(shù)值誤差定義如下

其中,?e表示精確解或細(xì)網(wǎng)格上得到的參考解.本文所有數(shù)值實(shí)驗(yàn)均在Visual Studio 2017+Intel Visual Fortran 的軟件平臺(tái)和CPU Xecon E3-1230+內(nèi)存12 GB 的臺(tái)式電腦上完成.

3.1 精度測(cè)試

考慮一維線性HJ 方程

以及二維線性HJ 方程

以上方程均考慮周期邊界條件,并分別計(jì)算至t=5 和t=0.5.表1 和表2 分別給出了一維和二維情形下由WCNS 格式和WENO 格式所得的L1和L∞數(shù)值誤差、精度階和CPU 運(yùn)行時(shí)間.由表可知,相比WENO 格式,WCNS 格式在光滑區(qū)域能達(dá)到所設(shè)計(jì)的7 階精度,其在精度和收斂性上明顯優(yōu)于經(jīng)典的WENO 格式.圖1 給出了兩種格式的L∞數(shù)值誤差與CPU 時(shí)間關(guān)系曲線圖.由此可知,當(dāng)獲得相同L∞數(shù)值誤差時(shí),WCNS 格式所耗的CPU 時(shí)間比WENO 格式少,因此計(jì)算效率略高.

表1 一維情形時(shí)兩種格式的數(shù)值誤差、精度階和CPU 運(yùn)行時(shí)間Table 1 Numerical errors,convergence rates and CPU time obtained with two schemes for the 1D case

表2 二維情形時(shí)兩種格式的數(shù)值誤差、精度階和CPU 運(yùn)行時(shí)間Table 2 Numerical errors,convergence rates and CPU time obtained with two schemes for the 2D case

圖1 兩種格式的計(jì)算誤差與CPU 時(shí)間關(guān)系曲線圖Fig.1 CPU times against numerical computing errors obtained with two schemes

3.2 分辨率測(cè)試

考慮滿足周期邊界條件的一維HJ 方程: ?t+?x=0,x∈[?1,1].本例測(cè)試兩種初值條件下格式的分辨率.第1 種初值條件為

網(wǎng)格數(shù)設(shè)置為N=100.采用7 階WCNS 格式和7 階WENO 格式計(jì)算本例.圖2 給出了基于初值條件(33)在t=11 時(shí)刻的數(shù)值結(jié)果,圖3 分別給出了基于初值條件(34) 在t=2 時(shí)刻和t=8 時(shí)刻的數(shù)值結(jié)果.由圖可知,WCNS 格式相比WENO 格式在間斷點(diǎn)附近具有更高的分辨率.在后面算例中僅采用7 階WCNS 格式計(jì)算.

圖2 基于初值條件(33)在 t=11 時(shí)刻所得數(shù)值解比較Fig.2 Comparisons of numerical solutions at time t=11 based on the initial condition (33)

圖3 基于初值條件(34)在不同時(shí)刻所得數(shù)值解比較Fig.3 Comparisons of numerical solutions at different times based on the initial condition (34)

3.3 一維凸和非凸Hamilton 問(wèn)題[11]

考慮一維具有凸Hamilton 函數(shù)的HJ 方程

以及一維具有非凸Hamilton 函數(shù)的HJ 方程

考慮周期邊界條件,并采用三種粗細(xì)不同的網(wǎng)格,即網(wǎng)格數(shù)設(shè)置為N=50,100,200.采用7階WCNS格式對(duì)以上兩個(gè)方程分別計(jì)算至此時(shí)兩個(gè)方程的解均會(huì)出現(xiàn)間斷.圖4 給出了基于兩個(gè)方程在不同網(wǎng)格下所得數(shù)值解.由圖可知,在3 種網(wǎng)格下所得數(shù)值解與參考解均吻合得很好.此外,WCNS 格式具有很好的間斷捕捉能力.

圖4 一維(a)凸和(b)非凸Hamilton 問(wèn)題的數(shù)值解Fig.4 Numerical solutions of 1D (a) convex and (b) nonconvex Hamilton problems

3.4 二維凸和非凸Hamilton 問(wèn)題

考慮二維具有凸Hamilton 函數(shù)的HJ 方程

和二維具有非凸Hamilton 函數(shù)的HJ 方程

考慮周期邊界條件,并采用3 種粗細(xì)不同的網(wǎng)格,即網(wǎng)格數(shù)為M×N=50×50,100×100,200×200.兩個(gè)方程分別計(jì)算至以測(cè)試WCNS 格式對(duì)間斷解的捕捉能力.圖5 和圖6 分別給出了基于兩個(gè)方程在50×50 網(wǎng)格下所得數(shù)值解的曲面圖和在各種網(wǎng)格下沿直線x=y的截面圖.由圖可知,WCNS格式對(duì)于該算例也具有很好的間斷捕捉能力.此外,在3 種網(wǎng)格下所得數(shù)值解與參考解均吻合得很好,因此,在后面算例中均考慮粗網(wǎng)格,即網(wǎng)格數(shù)為 5 0×50.

圖5 二維凸Hamilton 問(wèn)題的數(shù)值解(a)曲面圖和(b)截面圖Fig.5 (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D convex Hamilton problem

圖6 二維非凸Hamilton 問(wèn)題的數(shù)值解(a)曲面圖和(b)截面圖Fig.6 (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D nonconvex Hamilton problem

3.5 二維完全問(wèn)題[10]

求解周期邊界條件下二維HJ 方程

該問(wèn)題具有精確解

其中,x=q?tsinr,y=r+tcosq.當(dāng)t<1.0 時(shí),方程具有光滑解;當(dāng)t≥1.0 時(shí),方程的解會(huì)出現(xiàn)間斷.圖7給出了t=0.8 時(shí)刻和t=1.5 時(shí)刻數(shù)值解的曲面圖和沿直線x=y的截面圖.由圖可知,數(shù)值解和精確解吻合,在間斷附近WCNS 格式基本無(wú)振蕩.

圖7 二維完全問(wèn)題的數(shù)值解(a),(b) 曲面圖和(c),(d)截面圖Fig.7 (a),(b) Surfaces and (c),(d) cross-sectional diagrams of numerical solutions of the fully problem

3.6 二維最優(yōu)控制問(wèn)題 [13,15]

圖8 二維最優(yōu)控制問(wèn)題的數(shù)值解(a),(b)曲面圖和(c),(d)截面圖Fig.8 (a),(b) Surfaces and (c),(d) cross-sectional diagrams of the numerical solutions of the optimal control problem

圖8 二維最優(yōu)控制問(wèn)題的數(shù)值解(a),(b)曲面圖和(c),(d)截面圖 (續(xù))Fig.8 (a),(b)Surfaces and (c),(d) cross-sectional diagrams of the numerical solutions of the optimal control problem (continued)

求解周期邊界條件下二維最優(yōu)控制問(wèn)題的截面圖.由圖可知,WCNS 格式精度高、分辨率好.

3.7 二維傳播面問(wèn)題[19,29]

求解周期邊界條件下二維傳播面問(wèn)題

圖9 二維傳播面問(wèn)題的 ε=0 時(shí)數(shù)值解Fig.9 Numerical solutions of the propagating surface problem with ε=0

其中,ε >0 為一個(gè)很小的常數(shù),K為平均曲率了 ε=0 和 ε=0.1 時(shí)不同時(shí)刻下的數(shù)值解曲面圖.由圖可知,WCNS 格式在數(shù)值模擬該問(wèn)題時(shí)基本無(wú)振蕩,具有很好的分辨率.

圖10 二維傳播面問(wèn)題的 ε=0.1 時(shí)數(shù)值解Fig.10 Numerical solutions of the propagating surface problem with ε=0.1

4 結(jié)論

本文針對(duì)非穩(wěn)態(tài)HJ 方程設(shè)計(jì)了7 階精度WCNS 格式.采用單調(diào)的Lax-Friedrichs 型通量分裂方法計(jì)算Hamilton 數(shù)值通量,8 階精度的混合型中心差分格式近似節(jié)點(diǎn)處一階空間導(dǎo)數(shù)的左、右極限值,并通過(guò)基于七點(diǎn)模板的WENO 型非線性插值方法得到半節(jié)點(diǎn)函數(shù)值的高階逼近.3 階顯式TVD Runge-Kutta 時(shí)間離散方法用于求解空間離散化后的半離散格式.精度測(cè)試算例結(jié)果表明,本文提出的WCNS 格式能夠很好的模擬HJ 方程的精確解并且在光滑區(qū)域能夠達(dá)到設(shè)計(jì)的7 階精度.相比經(jīng)典的7 階WENO 格式,WCNS 格式在精度和收斂性方面更優(yōu),且獲得相同誤差時(shí),所耗CPU 時(shí)間更短,因此計(jì)算效率略高.分辨率測(cè)試算例結(jié)果表明,相比經(jīng)典的7 階WENO 格式,WCNS 格式具有更好的分辨率和更強(qiáng)的間斷捕捉能力.其他一維和二維測(cè)試算例結(jié)果均表明WCNS 格式精度高、分辨率高、間斷捕捉能力強(qiáng).下一步工作,進(jìn)一步優(yōu)化全局模板和子模板的光滑度量指標(biāo),提高WCNS 格式計(jì)算效率.