胡迎港 蔣艷群 黃曉倩

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

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

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

1 數(shù)值方法

1.1 一維HJ 方程的WCNS 格式

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

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

對式(10)右端在半節(jié)點xi+1/2處進行Taylor 展開,得到

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

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

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

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

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

1.2 二維HJ 方程的WCNS 格式

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

這里同樣考慮均勻網(wǎng)格,即網(wǎng)格節(jié)點設(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)改寫為如下形式

將其代入式(5)得到

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

3 數(shù)值實驗

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

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

3.1 精度測試

考慮一維線性HJ 方程

以及二維線性HJ 方程

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

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

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

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

3.2 分辨率測試

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

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

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

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

3.3 一維凸和非凸Hamilton 問題[11]

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

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

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

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

3.4 二維凸和非凸Hamilton 問題

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

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

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

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

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

3.5 二維完全問題[10]

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

該問題具有精確解

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

圖7 二維完全問題的數(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)控制問題 [13,15]

圖8 二維最優(yōu)控制問題的數(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)控制問題的數(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)控制問題的截面圖.由圖可知,WCNS 格式精度高、分辨率好.

3.7 二維傳播面問題[19,29]

求解周期邊界條件下二維傳播面問題

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

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

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

4 結(jié)論

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