HIGH ORDER WEIGHT COMPACR NONLINEAR SCHEME FOR TIME-DEPENDENT HAMILTON-JACOBI EQUATIONS
-
摘要: Hamilton-Jacobi (HJ) 方程是一类重要的非线性偏微分方程, 在物理学、流体力学、图像处理、微分几何、金融数学、最优化控制理论等方面有着广泛的应用. 由于HJ方程的弱解存在但不唯一, 且解的导数可能出现间断, 导致其数值求解具有一定的难度. 本文提出了非稳态HJ方程的7阶精度加权紧致非线性格式 (WCNS). 该格式结合了Hamilton函数的Lax-Friedrichs型通量分裂方法和一阶空间导数左、右极限值的高阶精度混合节点和半节点型中心差分格式. 基于7点全局模板和4个4点子模板推导了半节点函数值的高阶线性逼近和4个低阶线性逼近, 以及全局模板和子模板的光滑度量指标. 为避免间断附近数值解产生非物理振荡以及提高格式稳定性, 采用WENO型非线性插值方法计算半节点函数值. 时间离散采用3阶TVD型Runge-Kutta方法. 通过理论分析验证了WCNS格式对于光滑解具有最佳的7阶精度. 为方便比较, 经典的7阶WENO格式也被推广用于求解HJ方程. 数值结果表明, 本文提出的WCNS格式能够很好地模拟HJ方程的精确解, 且在光滑区域能够达到7阶精度; 与经典的同阶WENO格式相比, WCNS格式在精度、收敛性和分辨率方面更优, 计算效率略高.
-
关键词:
- Hamilton-Jacobi方程 /
- 7阶WCNS /
- WENO插值 /
- 高分辨率 /
- 间断捕捉能力
Abstract: Hamilton-Jacobi (HJ) equations are an important class of nonlinear partial differential equations. They are often used in various applications, such as physics, fluid mechanics, image processing, differential geometry, financial mathematics, optimal control theory, and so on. Because the weak solutions of the HJ equations exist but are not unique, and the spatial derivatives of the solutions may be discontinuous, numerical difficulties arise in numerical solutions of these equations. This paper presents a seventh-order weighted compact nonlinear scheme (WCNS) for the time-dependent HJ equations. This scheme is composed of the monotone Lax-Friedrichs flux splitting method for the Hamilton functions and the high-order hybrid cell-node and cell-edge central differencing for the left and right limits of first-order spatial derivatives in the numerical Hamilton functions. A high-order linear approximation scheme and four low-order linear approximation schemes for the unknowns at half nodes are derived based on a seven-point global stencil and four four-point sub-stencils, respectively. The smoothness indicators of the global stencil and four sub-stencils are also derived. In order to avoid non-physical oscillations of numerical solutions near the discontinuities and improve the numerical stability of the designed scheme, the WENO-type nonlinear interpolation technique is adopted to compute the unknowns at half nodes. The third-order TVD Runge-Kutta method is used for time discretization. The presented WCNS scheme is verified to have the optimal seventh order of accuracy for smooth solutions by theoretical analysis. For the sake of comparison, the classical seventh-order WENO scheme for solving hyperbolic conservation laws is also extended to solve the HJ equations. Numerical results show that the presented WCNS scheme can well simulate the exact solutions and can achieve seventh-order accuracy in smooth regions. Compared with the classical WENO scheme of the same order, the presented WCNS scheme has better accuracy, convergence and resolution, and its computational efficiency is slightly higher. -
引 言
$ d $ 维空间下非稳态Hamilton-Jacobi (HJ)方程定义如下$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + H\left( {{{\boldsymbol{x}}},t,\phi ,\nabla \phi } \right) = 0} \\ {\phi \left( {{{\boldsymbol{x}}},0} \right) = {\phi _0}\left( {{\boldsymbol{x}}} \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \end{array}} \right\} $$ (1) 其中,
$ t \gt 0 $ ,$ {{\boldsymbol{x}}} = {({x_1}, {x_2},\cdots {x_d})^{\rm{T}}} $ ,$ \nabla \phi {\text{ = (}}{\phi _{{x_1}}},{\phi _{{x_2}}},\cdots{\text{,}}{\phi _{{x_d}}}{{\text{)}}^{\rm{T}}} $ . 这类方程在物理学、流体力学、图像处理、微分几何、金融数学、最优化控制理论等诸多领域中有着广泛应用[1-4]. 其解的性质很特殊, 如弱解存在但不唯一, 即使初值和Hamilton函数充分光滑, 方程解的导数也可能出现间断[5]. Crandall等[6]首次提出了HJ方程的黏性解的概念, 并证明了黏性解的存在唯一性. HJ方程的性质与双曲守恒律方程十分相似, 因此, 可以将双曲守恒律方程的一些成熟数值方法[7-9]推广用于求解HJ方程.目前已有较多的高精度算法被推广应用于HJ方程. 如Osher等[10]通过选择最光滑的候选模板重构数值通量方法设计了HJ方程的二阶本质无振荡 (ENO) 格式. Jiang等[11]通过将所有候选模板的重构通量进行了非线性加权得到了HJ方程的5阶加权本质无振荡 (WENO) 格式. 相比于ENO格式, WENO格式在光滑区域能达到更高阶精度, 同时在间断区域恢复为ENO格式. Bryson等[12]设计了HJ方程的5阶中心加权本质无振荡(CWENO)格式. Qiu等[13]提出了HJ方程的基于Hermite多项式的5阶HWENO格式, 该格式在重构通量时具有更好的紧致性. Ha等[14]为改善格式的间断捕捉能力, 采用拉格朗日型指数多项式建立了一种新的6阶WENO格式. 相比于其他类型的WENO格式, 该格式具有更低的计算成本. Cheng等[15]基于经典的5阶WENO格式, 采用4个二次多项式的非线性凸组合重构数值通量, 提高了格式的精度 (在光滑区域能达到6阶) 和健壮性. Zhu等[16]基于原始的5阶HWENO格式, 通过采用大小不同的候选模板构造了HJ方程的更加简单有效的HWENO格式. Rathan[17]提出了一种新的5阶WENO格式求解HJ方程, 该格式通过设计L1范数型的非线性权重来提高精度和分辨率. Abedian[18]采用对称的5阶WENO格式求解HJ方程. Kim等[19]基于指数多项式构造了一种新的3阶WENO格式. 其特点是可以很好地分辨出光滑区域和间断区域.
高阶精度加权紧致非线性格式 (WCNS) 是另一种求解双曲守恒律方程的有效方法. Deng等[20]基于WENO格式的构造思想, 在紧致非线性格式[21] (CNS) 中引入非线性WENO插值技术, 提出了高阶精度WCNS格式. 该格式被证实具有高分辨率和强间断捕捉能力等良好特性. Nonomura等[22]和Zhang等[23]分别对WCNS格式进行了改进, 提出7阶和9阶精度的WCNS格式. 并通过数值试验证明了格式的高分辨率和强间断捕捉能力. 为了进一步提高WCNS格式在复杂流体计算中的应用, Deng等[24]构造了混合半节点和节点的加权紧致非线性格式 (HWCNS), 使得格式在间断区域有了更高的分辨率. Nonomura等[25]设计一种更加健壮的混合型WCNS格式. Wong等[26]提出局部耗散加权紧致格式 (WCNS-LD). 该格式在光滑区域使用耗散较小的中心插值模板, 在包含间断区域使用耗散较大的迎风插值模板以抑制非物理振荡. Zhao等[27]提出一种新的可调参数的光滑度量指标, 进一步提高了WCNS格式的间断捕捉能力和健壮性. 洪正等[28]采用5阶WCNS格式对各向同性湍流通过正激波的情形进行直接数值模拟. Jiang等[29]为HJ方程设计了5阶精度的WCNS格式, 该格式相比于同阶WENO格式具有更好的收敛性和分辨率. Hiejima[30]基于TENO插值思想, 建立WCNS-T格式. 相比传统的WCNS格式, WCNS-T格式在强不连续和高频间断的捕捉上更具优势. Jiang等[31]针对等熵Navier-Stokes方程组提出半隐式时间推进的WCNS格式, 避免显式WCNS格式严格的CFL稳定性限制. Ma等[32]构造了非线性紧致插值的WCNS格式, 数值验证新的WCNS格式可用于大涡模拟.
本文在Jiang等[29]提出的5阶精度WCNS格式基础上, 进一步构造了非稳态HJ方程的7阶精度WCNS格式. HJ方程的Hamilton数值通量采用具有单调性的Lax-Friedrichs方法计算, 一阶空间导数的左、右极限值采用高阶精度混合节点和半节点型的8阶中心差分格式计算, 半节点函数值采用7阶WENO型非线性插值方法计算. 3阶TVD Runge-Kutta方法[33]用于HJ方程的时间离散. 最后通过一系列一维和二维数值算例验证WCNS格式在精度, 分辨率和间断捕捉能力等方面的特性.
1. 数值方法
1.1 一维HJ方程的WCNS格式
考虑如下形式的一维非稳态HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + H\left( {{\phi _x}} \right) = 0} \\ {\phi \left( {x,0} \right) = {\phi _0}\left( x \right)} \\ \qquad x \in [a,b] \end{array} } \right\} $$ (2) 为简单起见, 采用均匀网格, 即网格节点设置为
${x_i} = a + i\Delta x\;(i = 0,1,\cdots,N)$ . 其中,$\Delta x{\text{ = (}}b - a{\text{)/}}N$ 为网格尺寸,$N$ 为网格数. 方程(2)的半离散格式为$$ \frac{{{\rm{d}}{\phi _i}(t)}}{{{\rm{d}}t}} = - \hat H(\phi _{x,i}^ + ,\phi _{x,i}^ - ) $$ (3) 其中,
${\phi _i}(t)$ 为节点${x_i}$ 处$\phi ({x_i},t)$ 的数值近似,$\phi _{x,i}^ - $ 和$\phi _{x,i}^ + $ 分别为一阶空间导数${\phi _x}$ 在节点${x_i}$ 处的左、右极限值,$\hat H$ 为Hamilton函数$H$ 的数值通量. 为了提高格式的稳定性, 采用具有单调性的全局Lax-Friedrichs型通量分裂方法计算$\hat H$ , 即$$ \hat H(\phi _{x,i}^ + ,\phi _{x,i}^ - ) = H\left(\frac{{\phi _{x,i}^ + + \phi _{x,i}^ - }}{2}\right) - \frac{\alpha }{2}(\phi _{x,i}^ + - \phi _{x,i}^ - )$$ (4) 其中,
$\alpha = {\max _{{\phi _x}}}\left| {{H_{{\phi _x}}}} \right|$ .采用混合节点和半节点型的8阶中心差分格式计算
${\phi _x}$ 在节点${x_i}$ 处的左、右极限值$\phi _{x,i}^ \mp $ , 即$$ \begin{split} \hat \phi _{x,i}^ \mp = &\dfrac{{265}}{{175\Delta x}}(\hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R} - \hat \phi _{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R}) - \frac{1}{{4\Delta x}}({\phi _{i + 1}} - {\phi _{i - 1}})+ \\ & \frac{1}{{100\Delta x}}({\phi _{i + 2}} - {\phi _{i - 2}}) - \frac{1}{{2100\Delta x}}({\phi _{i + 3}} - {\phi _{i - 3}}) \end{split} $$ (5) 其中,
$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R} $ 为半节点${x_{i{\text{ + }}1/2}}$ 处函数$ \phi $ 的左、右状态值, 由7阶WENO型非线性插值方法计算得到. 由于$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^L $ 和$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^R $ 关于$x = {x_i}$ 镜像对称, 因此, 本文只给出$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^L $ 的计算过程, 且为了推导简洁省略上标“L”. 7阶WENO型非线性插值过程描述如下.选择7点模板
$ S = \left\{ {{x_{i - 3}},{x_{i - 2}}, \cdots ,{x_{i + 3}}} \right\} $ . 基于该模板得到如下6次Lagrange插值多项式$$\hat{\phi}(x)=\sum_{l=i-3}^{i+3} \prod_{\substack{j=i-3 \\ j \neq l}}^{i+3} \frac{x-x_j}{x_l-x_j} \phi_l $$ (6) 由此可得半节点函数值
${\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}$ 的一个7阶线性逼近$$ \begin{split} {{\hat \phi }_{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} =& \frac{1}{{1024}}( - 5{\phi _{i - 3}} + 42{\phi _{i - 2}} - 175{\phi _{i - 1}}+ \\ & 700{\phi _i} + 525{\phi _{i + 1}} - 70{\phi _{i + 2}} + 7{\phi _{i + 3}}) \end{split} $$ (7) 对式(7)右端在半节点
${x_{i{\text{ + }}1/2}}$ 处进行Taylor展开, 得到$$ {\hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = {\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} + A\Delta {x^7} + O(\Delta {x^8})$$ (8) 其中,
$A$ 为与$\Delta x$ 无关的常系数.全局模板
$ S $ 可划分为4个4点子模板, 即$ {S_k} = \left\{ {{x_{i + k - 3}},{x_{i + k - 2}},{x_{i + k - 1}},{x_{i + k}}} \right\} $ ,$ k = 0,1,2,3 $ . 在每个子模板$ {S_k} $ 上可推导出一个3次Lagrange插值多项式$$ \hat{\phi}_k(x)=\sum_{l=i+k-3}^{i+k} \prod_{\substack{j=i+k-3 \\ j \neq l}}^{i+k} \frac{x-x_j}{x_l-x_j} \phi_l,\;\;\;\; k=0,1,2,3$$ (9) 由此可得半节点函数值
${\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}$ 的4个4阶线性逼近$$ \left.\begin{array}{l} {{\hat \phi }_{0,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \dfrac{1}{{16}}( - 5{\phi _{i - 3}} + 21{\phi _{i - 2}} - 35{\phi _{i - 1}} + 35{\phi _i}) \\ {{\hat \phi }_{1,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \dfrac{1}{{16}}({\phi _{i - 2}} - 5{\phi _{i - 1}} + 15{\phi _i} + 5{\phi _{i + 1}}) \\ {{\hat \phi }_{2,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \dfrac{1}{{16}}( - {\phi _{i - 1}} + 9{\phi _i} + 9{\phi _{i + 1}} - {\phi _{i + 2}}) \\ {{\hat \phi }_{3,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \dfrac{1}{{16}}(5{\phi _i} + 15{\phi _{i + 1}} - 5{\phi _{i + 2}} + {\phi _{i + 3}}) \end{array}\right\} $$ (10) 对式(10)右端在半节点
${x_{i{\text{ + }}1/2}}$ 处进行Taylor展开, 得到$$ {\hat \phi _{k,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = {\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} + {B_k}\Delta {x^4} + O(\Delta {x^5}) $$ (11) 其中,
${B_k}(k = 0,1,2,3)$ 为与$\Delta x$ 无关的常系数. 半节点函数值$ {\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} $ 的高阶线性逼近式(7)可由式(10)中4个低阶线性逼近的线性组合得到$$ {\hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \sum\limits_{k = 0}^3 {{d_k}{{\hat \phi }_{k,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}} $$ (12) 其中,
${d_0} = {1 \mathord{\left/ {\vphantom {1 {64}}} \right. } {64}}$ ,${d_1} = {{21} \mathord{\left/ {\vphantom {{21} {64}}} \right. } {64}}$ ,${d_2} = {{35} \mathord{\left/ {\vphantom {{35} {64}}} \right. } {64}}$ ,${d_3} = {7 \mathord{\left/ {\vphantom {7 {64}}} \right. } {64}}$ 为理想线性权重.若直接采用线性格式(12)计算
$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{} $ , 则在间断附近数值解会出现非物理振荡. 为避免非物理振荡以及提高格式的数值稳定性, 基于WENO插值思想[7-9], 采用非线性权重$ {\omega _k} $ 代替线性权重${d_k}$ , 得到$$ {\hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \sum\limits_{k = 0}^3 {{\omega _k}{{\hat \phi }_{k,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}} $$ (13) 当所有子模板都处于光滑区域时, 非线性权重等于线性权重, 得到7阶精度的半节点函数逼近; 当某些子模板包含间断区域时, 令其非线性权重趋于零, 防止插值穿过间断区域, 最后得到4阶精度的半节点函数逼近. 本文选择WENO-Z型非线性权重, 定义如下
$$ {\omega _k} = \frac{{{\alpha _k}}}{{\displaystyle\sum\limits_{k = 0}^3 {{\alpha _k}} }},\;\;\;\;{\alpha _k} = {d_k}\left(1 + \frac{\tau }{{I{S_k} + \varepsilon }}\right) $$ (14) 式中,
$k = 0,1,2,3$ , 参数$\varepsilon $ 取值为小量${10^{ - 20}}$ , 以避免分母为0. 为了使得所设计WCNS格式在光滑区域能达到7阶精度, 全局光滑度量指标$\tau $ 定义为$ \tau = \left| {I{S_3} - I{S_2} - I{S_1} + I{S_0}} \right| $ . 其中,$I{S_k}$ 是子模板${S_k}$ $(k = 0,1,2, 3)$ 的光滑度量指标, 定义如下$$ I{S_k} = \sum\limits_{m = 1}^3 {{{[{{(\Delta x)}^m}\hat \phi _k^{(m)}({x_i})]}^2}} $$ (15) 其中,
$ \hat \phi _k^{(m)} $ 表示子模板${S_k}$ 上Language型插值多项式的$m$ 阶导数. 各个子模板的光滑度量指标$I{S_k}$ 的具体表达式如下$$ \begin{split} I{S_0} =& {\left( - \frac{1}{3}{\phi _{i - 3}} + \frac{3}{2}{\phi _{i - 2}} - 3{\phi _{i - 1}} + \frac{{11}}{6}\phi \right)^2}+ \\ & {( - {\phi _{i - 3}} + 4{\phi _{i - 2}} - 5{\phi _{i - 1}} + 2\phi )^2} +\\ & {( - {\phi _{i - 3}} + 3{\phi _{i - 2}} - 3{\phi _{i - 1}} + {\phi _i})^2}\\ \\ I{S_1} =& {\left(\frac{1}{6}{\phi _{i - 2}} - {\phi _{i - 1}} + \frac{1}{2}{\phi _i} + \frac{1}{3}{\phi _{i + 1}}\right)^2} +\\ & {({\phi _{i - 1}} - 2{\phi _i} + {\phi _{i + 1}})^2} + \\ & {( - {\phi _{i - 2}} + 3{\phi _{i - 1}} - 3{\phi _i} + {\phi _{i + 1}})^2}\\ \\ I{S_2} = &{\left( - \frac{1}{3}{\phi _{i - 1}} + 9{\phi _i} + 9{\phi _{i + 1}} - {\phi _{i + 2}}\right)^2} +\\ & {({\phi _{i - 1}} - 2{\phi _i} + {\phi _{i + 1}})^2} +\\ & {( - {\phi _{i - 1}} + 3{\phi _i} - 3{\phi _{i + 1}} + {\phi _{i + 2}})^2} \\ I{S_3} = &{\left( - \frac{{11}}{6}{\phi _i} + 3{\phi _{i + 1}} - \frac{3}{2}{\phi _{i + 2}} + \frac{1}{3}{\phi _{i + 3}}\right)^2} +\\ & {(2{\phi _i} - 5{\phi _{i + 1}} + 4{\phi _{i + 2}} - 3{\phi _{i + 3}})^2} +\\ & {( - {\phi _i} + 3{\phi _{i + 1}} - 3{\phi _{i + 2}} + {\phi _{i + 3}})^2} \end{split} $$ (16) 对于半离散格式(3), 采用如下3阶显式TVD Runge-Kutta 方法[33]求解
$$ \left.\begin{gathered} {\phi ^{\left( 1 \right)}} = {\phi ^n} - \Delta t\hat H(\phi _x^{n, + },\phi _x^{n, - }) \\ {\phi ^{\left( 2 \right)}} = \frac{3}{4}{\phi ^n} + \frac{1}{4}{\phi ^{\left( 1 \right)}} - \frac{1}{4}\Delta t\hat H(\phi _x^{\left( 1 \right), + },\phi _x^{\left( 1 \right), - }) \\ {\phi ^{(n + 1)}} = \frac{1}{3}{\phi ^n} + \frac{2}{3}{\phi ^{\left( 1 \right)}} - \frac{2}{3}\Delta t\hat H(\phi _x^{\left( 2 \right), + },\phi _x^{\left( 2 \right), - }) \end{gathered}\right\} $$ (17) 其中, 时间步长
$\Delta t = {t^{n + 1}} - {t^n}$ 需满足CFL稳定性条件$$ \Delta t \leqslant C\frac{{\Delta x}}{\alpha } $$ (18) 其中,
$C$ 为CFL数, 本文取值为0.6.本文将文献[9]中求解双曲守恒律方程的7阶WENO格式推广用于求解HJ方程. 式(4)中一阶空间导数
${\phi _x}$ 在节点${x_i}$ 处的左、右极限值$\phi _{x,i}^ \mp $ 的计算格式如下$$ \hat \phi _{x,i}^ \mp = \frac{1}{{\Delta x}}(\hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R} - \hat \phi _{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R}) $$ (19) 半节点
${x_{i{\text{ + }}1/2}}$ 处函数$ \phi $ 的左、右状态值$ \hat \phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}^{L,R} $ 由基于7点模板$ S = \left\{ {{x_{i - 3}},{x_{i - 2}}, \cdots ,{x_{i + 3}}} \right\} $ 的7阶WENO重构方法[9]计算得到. 由此可见, 本文设计的7阶WCNS格式与经典的7阶WENO格式采用相同的全局模板$ S $ .1.2 二维HJ方程的WCNS格式
考虑如下形式的二维非稳态HJ方程
$$ \left. \begin{array}{*{20}{c}} {{\phi _t} + H\left( {{\phi _x},{\phi _y}} \right) = 0} \\ {\phi \left( {x,y,0} \right) = {\phi _0}\left( {x,y} \right)} \\ \qquad (x,y) \in [a,b] \times [c,d] \end{array} \right\} $$ (20) 这里同样考虑均匀网格, 即网格节点设置为
${x_i} = i\Delta x\;(i = 0,1,\cdots,M)$ ,${y_j} = j\Delta y\;(j = 0,1,\cdots,N)$ . 其中,$x$ 和$y$ 方向的网格尺寸分别为$\Delta x{\text{ = (}}b - a{\text{)/}}M$ ,$\Delta y{\text{ = (}}d - c{\text{)/}}N$ , 网格数为$M \times N$ .方程(20)的半离散格式为
$$ \frac{{{\rm{d}}{\phi _{ij}}(t)}}{{{\rm{d}}t}} = - \hat H(\phi _{x,ij}^ + ,\phi _{x,ij}^ - ,\phi _{y,ij}^ + ,\phi _{y,ij}^ - ) $$ (21) 其中,
${\phi _{ij}}(t)$ 为函数$\phi ({x_i},{y_j},t)$ 在网格节点$({x_i},{y_j})$ 处的数值近似,$ \phi _{x,ij}^ \mp $ 和$ \phi _{y,ij}^ \mp $ 分别为一阶偏导数${\phi _x}$ 和${\phi _y}$ 在$({x_i},{y_j})$ 点处的左、右极限值. Hamilton 函数$H$ 的Lax-Friedrichs型数值通量为$$ \begin{split}& \hat H(\phi _{x,ij}^ + ,\phi _{x,ij}^ - ,\phi _{y,ij}^ + ,\phi _{y,ij}^ - ) = H\left(\frac{{\phi _{x,ij}^ + + \phi _{x,ij}^ - }}{2},\frac{{\phi _{y,ij}^ + + \phi _{y,ij}^ - }}{2}\right)- \\ &\qquad \frac{{{\alpha _x}}}{2}(\phi _{x,ij}^ + - \phi _{x,ij}^ - ) - \frac{{{\alpha _y}}}{2}(\phi _{y,ij}^ + - \phi _{y,ij}^ - ) \end{split} $$ (22) 其中,
${\alpha _x} = {\max _{{\phi _x}}}\left| {{H_{{\phi _x}}}} \right|,{\alpha _y} = {\max _{{\phi _y}}}\left| {{H_{{\phi _y}}}} \right|$ . 采用逐维形式的WCNS格式在$x$ 和$y$ 方向上分别求解偏导数${\phi _x}$ 和${\phi _y}$ 的左、右极限值$\phi _{x,ij}^ \mp $ 和$\phi _{y,ij}^ \mp $ . 半离散格式(21)同样采用3阶显式TVD Runge-Kutta 方法计算, 时间步长$\Delta t$ 需满足$$ \Delta t \leqslant C\frac{1}{{\max ({\alpha _x}/\Delta x,{\alpha _y}/\Delta y)}} $$ (23) 2. WCNS格式的精度分析
基于一维HJ方程, 式(13)改写为如下形式
$$ \begin{split} & {{\hat \phi }_{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} = \sum\limits_{k = 0}^3 {{d_k}{{\hat \phi }_{k,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} + } \sum\limits_{k = 0}^3 {({\omega _k} - {d_k}){{\hat \phi }_{k,i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}} = \\ &\qquad {\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} + A\Delta {x^7} + O(\Delta {x^8}) +\\ &\qquad \displaystyle\sum\limits_{k = 0}^3 {({\omega _k} - {d_k})({\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} + {B_k}\Delta {x^4} + O(\Delta {x^5}))} \end{split} $$ (24) 将其代入式(5)得到
$$ \begin{split} {{\hat \phi }_{x,i}} = &{\phi _{x,i}} + O(\Delta {x^8})+ \\ & \frac{{265}}{{175}}\Biggr\{\frac{1}{{\Delta x}}\Biggr[\sum\limits_{k = 0}^3 {(\omega _k^ + - {d_k}){\phi _{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}} - \sum\limits_{k = 0}^3 {(\omega _k^ - - {d_k}){\phi _{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}} } \Biggr]+ \\ & \Delta {x^3}\sum\limits_{k = 0}^3 {{B_k}(\omega _k^ + - \omega _k^ - )} + \sum\limits_{k = 0}^3 {(\omega _k^ + - \omega _k^ - )O(\Delta {x^4})}\Biggr\} \end{split} $$ (25) 因此, WCNS格式达到最佳7阶精度的充分条件为
$$\left. \begin{array}{l} \displaystyle\sum\limits_{k = 0}^3 {(\omega _k^ \pm - {d_k})} \leqslant O(\Delta {x^8}) \\ \omega _k^ \pm - {d_k} \leqslant O(\Delta {x^4}) \end{array} \right\}$$ (26) 由于权重满足
$ \displaystyle\sum\limits_{k = 0}^3 {\omega _k^ \pm } = \displaystyle\sum\limits_{k = 0}^3 {{d_k} = 1} $ , 第一个条件总是满足的. 将子模板$ {S_k} $ 的光滑度量指标$I{S_k}$ 在$ {x_i} $ 点进行Taylor展开, 得到 (省略下标$i$ )$$ \left.\begin{gathered} I{S_0} = {({\phi ^{(1)}}\Delta x)^2} + {\phi ^{(2)}}^2\Delta {x^4} - \frac{1}{2}{\phi ^{(1)}}{\phi ^{(4)}}\Delta {x^5} +\\ \left( - \frac{{11}}{6}{\phi ^{(2)}}{\phi ^{(4)}} + {\phi ^{(3)}}^2\right)\Delta {x^6} + O(\Delta {x^7}) \\ I{S_1} = {({\phi ^{(1)}}\Delta x)^2} + {\phi ^{(2)}}^2\Delta {x^4} + \frac{1}{6}{\phi ^{(1)}}{\phi ^{(4)}}\Delta {x^5} +\\ \left(\frac{1}{6}{\phi ^{(2)}}{\phi ^{(4)}} + {\phi ^{(3)}}^2\right)\Delta {x^6} + O(\Delta {x^7})\\ I{S_2} = {({\phi ^{(1)}}\Delta x)^2} + {\phi ^{(2)}}^2\Delta {x^4} - \frac{1}{6}{\phi ^{(1)}}{\phi ^{(4)}}\Delta {x^5}+ \\ \left(\frac{1}{6}{\phi ^{(2)}}{\phi ^{(4)}} + {\phi ^{(3)}}^2\right)\Delta {x^6} + O(\Delta {x^7}) \\ I{S_3} = {({\phi ^{(1)}}\Delta x)^2} + {\phi ^{(2)}}^2\Delta {x^4} + \frac{1}{2}{\phi ^{(1)}}{\phi ^{(4)}}\Delta {x^5} +\\ \left( - \frac{{11}}{6}{\phi ^{(2)}}{\phi ^{(4)}} + {\phi ^{(3)}}^2\right)\Delta {x^6} + O(\Delta {x^7}) \end{gathered}\right\} $$ (27) 即
$$ I{S_k} = {({\phi ^{(1)}}\Delta x)^2}(1 + O(\Delta {x^2})) $$ (28) 代入式(14)得到
$$ \frac{{{\omega _k}}}{{{d_k}}} = 1 + \frac{\tau }{{I{S_k} + \varepsilon }} \leqslant 1 + O(\Delta {x^4})$$ (29) 于是, 当
$\tau \leqslant O(\Delta {x^6})$ 时, WCNS格式可以达到7阶精度. 由式(27)可知$ \tau = O(\Delta {x^6}) $ , 因此, 本文提出的WCNS格式能达到最佳7阶精度.3. 数值实验
本节数值测试7阶WCNS格式的性能. 在精度测试算例中,
$\Delta t = {(\Delta x)^{7/3}}$ ,$ {L_1} $ 和$ {L_\infty } $ 范数数值误差定义如下$$ \left.\begin{gathered} {L_1} = \sum\limits_{i = 0}^N {\frac{{\left| {{\phi _i} - {\phi _e}} \right|}}{{N + 1}}} \\ {L_\infty } = \mathop {\max }\limits_{0 \leqslant i \leqslant N} \left| {{\phi _i} - {\phi _e}} \right| \end{gathered} \right\}$$ (30) 其中,
${\phi _e}$ 表示精确解或细网格上得到的参考解. 本文所有数值实验均在Visual Studio 2017 + Intel Visual Fortran的软件平台和CPU Xecon E3-1230 + 内存12 GB 的台式电脑上完成.3.1 精度测试
考虑一维线性HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + {\phi _x} = 0} \\ {\phi (x,0) = \sin ({\text{π}} x)} \\ \qquad x \in [ - 1,1] \end{array}} \right\} $$ (31) 以及二维线性HJ方程
$$ \left. \begin{array}{*{20}{l}} {{\phi _t} + {\phi _x} + {\phi _y} = 0} \\ {\phi (x,y,0) = \sin [{\text{π}} (x + y)]} \\ \qquad (x,y)\in {{[ - 1,1]}^2} \end{array} \right\} $$ (32) 以上方程均考虑周期边界条件, 并分别计算至
$ t = 5 $ 和$ t = 0.5 $ . 表1和表2分别给出了一维和二维情形下由WCNS格式和WENO格式所得的${L_1}$ 和${L_\infty }$ 数值误差、精度阶和CPU运行时间. 由表可知, 相比WENO格式, WCNS格式在光滑区域能达到所设计的7阶精度, 其在精度和收敛性上明显优于经典的WENO格式. 图1给出了两种格式的${L_\infty }$ 数值误差与CPU时间关系曲线图. 由此可知, 当获得相同${L_\infty }$ 数值误差时, WCNS格式所耗的CPU时间比WENO格式少, 因此计算效率略高.表 1 一维情形时两种格式的数值误差、精度阶和CPU运行时间Table 1. Numerical errors, convergence rates and CPU time obtained with two schemes for the 1D caseMethod $N$ ${L_1}$error Order ${L_\infty }$error Order CPU time WCNS 60 2.00 × 10−8 — 3.72 × 10−8 — 0.078125 80 2.67 × 10−9 7.00 4.90 × 10−9 7.04 0.218750 100 5.60 × 10−10 7.00 1.02 × 10−9 7.05 0.453125 120 1.56 × 10−10 7.00 2.81 × 10−10 7.05 0.796875 140 5.33 × 10−11 6.98 9.52 × 10−11 7.03 1.312500 WENO 60 7.73 × 10−8 — 5.94 × 10−7 — 0.078125 80 1.33 × 10−8 6.13 1.31 × 10−7 5.24 0.203125 100 3.36 × 10−9 6.15 4.09 × 10−8 5.23 0.453125 120 1.14 × 10−9 5.94 1.57 × 10−8 5.24 0.828125 140 4.56 × 10−10 5.93 6.97× 10−9 5.27 1.093750 表 2 二维情形时两种格式的数值误差、精度阶和CPU运行时间Table 2. Numerical errors, convergence rates and CPU time obtained with two schemes for the 2D caseMethod $M \times N$ ${L_1}$error Order ${L_\infty }$error Order CPU time WCNS 60 × 60 2.95× 10−8 — 4.80× 10−8 — 1.578125 80 × 80 3.95× 10−9 6.99 6.42× 10−9 6.99 4.843750 100 × 100 8.28 × 10−10 7.00 1.34× 10−9 7.01 12.65625 120 × 120 2.31 × 10−10 7.00 3.74 × 10−10 7.02 27.85937 140 × 140 7.86 × 10−11 7.00 1.27 × 10−10 7.01 54.92187 WENO 60 × 60 2.96× 10−8 — 1.88 × 10−7 — 1.140625 80 × 80 4.48× 10−9 6.57 3.84× 10−8 5.52 3.765625 100 × 100 1.07× 10−9 6.41 1.11× 10−8 5.55 9.984375 120 × 120 3.33 × 10−10 6.41 3.94× 10−9 5.70 22.09375 140 × 140 1.24 × 10−10 6.41 1.80× 10−9 5.09 43.14062 ![]()
图 1 两种格式的计算误差与CPU时间关系曲线图Figure 1. CPU times against numerical computing errors obtained with two schemes3.2 分辨率测试
考虑满足周期边界条件的一维HJ方程:
$ {\phi _t} + {\phi _x} = 0, x \in [ - 1,1] $ . 本例测试两种初值条件下格式的分辨率. 第1种初值条件为$$ \phi \left( {x,0} \right) = \left\{ {\begin{array}{*{20}{l}} { - x\sin \left( {{{3{\text{π}} {x^2}} \mathord{\left/ {\vphantom {{3{\text{π}} {x^2}} 2}} \right. } 2}} \right),{\text{ }}x \in \left[ { - 1, - {1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}} \right)} \\ {\left| {\sin \left( {2{\text{π}} x} \right)} \right|,{\text{ }}x \in \left[ { - {1 \mathord{\left/ {\vphantom {1 3}} \right. } 3},{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}} \right)} \\ {2x - 1 - {{\sin \left( {3{\text{π}} x} \right)} \mathord{\left/ {\vphantom {{\sin \left( {3{\text{π}} x} \right)} 6}} \right. } 6},{\text{ }}x \in \left[ {{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3},1} \right]} \end{array}} \right. $$ (33) 第2种初值条件为:
$\phi \left( {x,0} \right) = g\left( {x - 0.5} \right)$ . 其中,$g\left( x \right)$ 的具体表达式如下$$ \begin{gathered} g\left( x \right) = - \left( {{{\sqrt 3 } \mathord{\left/ {\vphantom {{\sqrt 3 } 2}} \right. } 2} + {9 \mathord{\left/ {\vphantom {9 2}} \right. } 2} + {{2{\text{π}} } \mathord{\left/ {\vphantom {{2{\text{π}} } 3}} \right. } 3}} \right)\left( {x + 1} \right) + \\ \left\{ {\begin{array}{*{20}{l}} {2\cos \left( {{{3{\text{π}} {x^2}} \mathord{\left/ {\vphantom {{3{\text{π}} {x^2}} 2}} \right. } 2}} \right) - \sqrt 3 ,{\text{ }}x \in \left[ { - 1, - {1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}} \right)} \\ {{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2} + 3\cos \left( {2{\text{π}} x} \right),{\text{ }}x \in \left[ { - {1 \mathord{\left/ {\vphantom {1 3}} \right. } 3},0} \right)} \\ {{{15} \mathord{\left/ {\vphantom {{15} 2}} \right. } 2} - 3\cos \left( {2{\text{π}} x} \right),{\text{ }}x \in \left[ {0,{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}} \right)} \\ {{{\left[ {28 + 4{\text{π}} + \cos \left( {3{\text{π}} x} \right)} \right]} \mathord{\left/ {\vphantom {{\left( {28 + 4{\text{π}} + \cos \left( {3{\text{π}} x} \right)} \right)} 3}} \right. } 3} + 6{\text{π}} \left( {{x^2} - x} \right),x \in \left[ {{1 \mathord{\left/ {\vphantom {1 3}} \right. } 3},1} \right]} \end{array}} \right. \end{gathered} $$ (34) 网格数设置为
$N = 100$ . 采用7阶WCNS格式和7阶WENO格式计算本例. 图2给出了基于初值条件(33)在$t = 11$ 时刻的数值结果, 图3分别给出了基于初值条件(34)在$t = 2$ 时刻和$t = 8$ 时刻的数值结果. 由图可知, WCNS格式相比WENO格式在间断点附近具有更高的分辨率. 在后面算例中仅采用7阶WCNS格式计算.![]()
图 2 基于初值条件(33)在$t = 11$ 时刻所得数值解比较Figure 2. Comparisons of numerical solutions at time$t = 11$ based on the initial condition (33)![]()
图 3 基于初值条件(34)在不同时刻所得数值解比较Figure 3. Comparisons of numerical solutions at different times based on the initial condition (34)3.3 一维凸和非凸Hamilton问题[11]
考虑一维具有凸Hamilton函数的HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + \dfrac{{{{({\phi _x} + 1)}^2}}}{2} = 0} \\ {\phi (x,0) = - \cos ({\text{π}} x)} \\ \qquad x \in [ - 1,1] \end{array}} \right\}$$ (35) 以及一维具有非凸Hamilton函数的HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} - \cos ({\phi _x} + 1) = 0} \\ {\phi (x,0) = - \cos ({\text{π}} x)} \\ \qquad x \in [ - 1,1] \end{array} } \right\} $$ (36) 考虑周期边界条件, 并采用三种粗细不同的网格, 即网格数设置为
$N = 50,100,200$ . 采用7阶WCNS格式对以上两个方程分别计算至$t = {{3.5} \mathord{\left/ {\vphantom {{3.5} {{{\text{π}} ^2}}}} \right. } {{{\text{π}} ^2}}}$ 和$t = {{1.5} \mathord{\left/ {\vphantom {{1.5} {{{\text{π}} ^2}}}} \right. } {{{\text{π}} ^2}}}$ , 此时两个方程的解均会出现间断. 图4给出了基于两个方程在不同网格下所得数值解. 由图可知, 在3种网格下所得数值解与参考解均吻合得很好. 此外, WCNS格式具有很好的间断捕捉能力.![]()
图 4 一维(a)凸和(b)非凸Hamilton问题的数值解Figure 4. Numerical solutions of 1D (a) convex and (b) nonconvex Hamilton problems3.4 二维凸和非凸Hamilton问题
考虑二维具有凸Hamilton函数的HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + \dfrac{{{{({\phi _x} + {\phi _y} + 1)}^2}}}{2} = 0 } \\ {\phi (x,0) = - \cos \left[{\text{π}} \left(\dfrac{{x + y}}{2}\right)\right]}\\ \qquad (x,y) \in {[ - 1,1]^2} \end{array}} \right\}$$ (37) 和二维具有非凸Hamilton函数的HJ方程
$$ \left. \begin{array}{*{20}{l}} {{\phi _t} - \cos ({\phi _x} + {\phi _y} + 1) = 0} \\ {\phi (x,0) = - \cos \left[{\text{π}} \left(\dfrac{{x + y}}{2}\right)\right]}\\ \qquad (x,y) \in {{[ - 1,1]}^2} \end{array} \right\} $$ (38) 考虑周期边界条件, 并采用3种粗细不同的网格, 即网格数为
$M \times N = 50 \times 50,100 \times 100,200 \times 200$ . 两个方程分别计算至$t = {{3.5} \mathord{\left/ {\vphantom {{3.5} {{{\text{π}} ^2}}}} \right. } {{{\text{π}} ^2}}}$ 和$t = {{1.5} \mathord{\left/ {\vphantom {{1.5} {{{\text{π}} ^2}}}} \right. } {{{\text{π}} ^2}}}$ , 以测试WCNS格式对间断解的捕捉能力. 图5和图6分别给出了基于两个方程在50×50网格下所得数值解的曲面图和在各种网格下沿直线$x = y$ 的截面图. 由图可知, WCNS格式对于该算例也具有很好的间断捕捉能力. 此外, 在3种网格下所得数值解与参考解均吻合得很好, 因此, 在后面算例中均考虑粗网格, 即网格数为$50 \times 50$ .![]()
图 5 二维凸Hamilton问题的数值解(a)曲面图和(b)截面图Figure 5. (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D convex Hamilton problem![]()
图 6 二维非凸Hamilton问题的数值解(a)曲面图和(b)截面图Figure 6. (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D nonconvex Hamilton problem3.5 二维完全问题[10]
求解周期边界条件下二维HJ方程
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} + {\phi _x}{\phi _y} = 0 } \\ {{\phi _0}(x,y,0) = \sin x + \cos y}\\ \qquad (x,y) \in {[ - {\text{π}} ,{\text{π}} ]^2} \end{array}} \right\}$$ (39) 该问题具有精确解
$$ \phi \left( {x,y,t} \right) = - \cos q \sin r + \sin q + \cos r $$ 其中,
$x = q - t\sin r ,\;\;y = r + t\cos q $ . 当$t \lt 1.0$ 时, 方程具有光滑解; 当$t \geqslant 1.0$ 时, 方程的解会出现间断. 图7给出了$t = 0.8$ 时刻和$t = 1.5$ 时刻数值解的曲面图和沿直线$x = y$ 的截面图. 由图可知, 数值解和精确解吻合, 在间断附近WCNS格式基本无振荡.3.6 二维最优控制问题 [13,15]
求解周期边界条件下二维最优控制问题
$$ \left. \begin{array}{*{20}{l}} \Big[ {\phi _t} + {\phi _x}\sin y + (\sin x + {\rm{sign}}({\phi _y})){\phi _y} - \\ \qquad \dfrac{1}{2}{\sin ^2}y - 1 + \cos (x) \Big] = 0 \\ {\phi (x,y,0) = 0} \\\qquad (x,y) \in {{[ - {\text{π}} ,{\text{π}} ]}^2}\end{array} \right\} $$ (40) 图8给出了
$t = 1.0$ 时刻方程的数值解曲面图和控制方程$\omega = {\rm{sign}}({\phi _y})$ 的数值解曲面图 (注: 在网格点$({x_i},{y_j})$ 处${\phi _{y,ij}} = (\phi _{y,ij}^ - + \phi _{y,ij}^ + )/2$ ) , 以及两者沿直线$x = 2.5$ 的截面图. 由图可知, WCNS格式精度高、分辨率好.![]()
图 8 二维最优控制问题的数值解(a), (b)曲面图和(c), (d)截面图Figure 8. (a), (b) Surfaces and (c), (d) cross-sectional diagrams of the numerical solutions of the optimal control problem![]()
8 二维最优控制问题的数值解(a), (b)曲面图和(c), (d)截面图 (续)8. (a), (b)Surfaces and (c), (d) cross-sectional diagrams of the numerical solutions of the optimal control problem (continued)![]()
7 二维完全问题的数值解(a), (b) 曲面图和(c), (d)截面图7. (a), (b) Surfaces and (c), (d) cross-sectional diagrams of numerical solutions of the fully problem3.7 二维传播面问题[19,29]
求解周期边界条件下二维传播面问题
$$ \left. {\begin{array}{*{20}{l}} {{\phi _t} - (1 - \varepsilon K)\sqrt {\phi _x^2 + \phi _y^2 + 1} = 0} \\ {\phi (x,y,0) = 1 - 0.25[\cos (2{\text{π}} x) - 1][\cos (2{\text{π}} y) - 1]} \\ \qquad (x,y) \in {[0,1]^2} \end{array}} \right\} $$ (41) 其中,
$\varepsilon \gt 0$ 为一个很小的常数,$K$ 为平均曲率$$ K = - \frac{{{\phi _{xx}}(\phi _x^2 + 1) + {\phi _{yy}}(\phi _x^2 + 1) - {\phi _{xy}}{\phi _x}{\phi _y}}}{{{{(\phi _x^2 + \phi _y^2 + 1)}^{1.5}}}} $$ 式中, 一阶导数
${\phi _x}$ 和${\phi _y}$ 的计算方法为${\phi _{x,ij}} = (\phi _{x,ij}^ - + \phi _{x,ij}^ + )/2$ ,${\phi _{y,ij}} = (\phi _{y,ij}^ - + \phi _{y,ij}^ + )/2$ ; 2阶偏导数${\phi _{xx,ij}}$ 和${\phi _{yy,ij}}$ 分别由$x$ 和$y$ 方向上基于7点模板的插值多项式关于$x$ 和$y$ 二次求导得到; 混合2阶偏导数${\phi _{xy}}$ 采用逐维方法计算, 即先固定$y$ 方向, 由$x$ 方向上插值多项式求导得到${\phi _{x,ij}}$ , 然后在模板各个节点处的$y$ 方向上由插值多项式求导得到${\phi _{xy,ij}}$ . 图9和图10分别给出了$\varepsilon = 0$ 和$\varepsilon = 0.1$ 时不同时刻下的数值解曲面图. 由图可知, WCNS格式在数值模拟该问题时基本无振荡, 具有很好的分辨率.![]()
图 9 二维传播面问题的$\varepsilon = 0$ 时数值解Figure 9. Numerical solutions of the propagating surface problem with$\varepsilon = $ 0![]()
图 10 二维传播面问题的$\varepsilon = 0.1$ 时数值解Figure 10. Numerical solutions of the propagating surface problem with$\varepsilon = $ 0.14. 结 论
本文针对非稳态HJ方程设计了7阶精度WCNS格式. 采用单调的Lax-Friedrichs型通量分裂方法计算Hamilton数值通量, 8阶精度的混合型中心差分格式近似节点处一阶空间导数的左、右极限值, 并通过基于七点模板的WENO型非线性插值方法得到半节点函数值的高阶逼近. 3阶显式TVD Runge-Kutta时间离散方法用于求解空间离散化后的半离散格式. 精度测试算例结果表明, 本文提出的WCNS格式能够很好的模拟HJ方程的精确解并且在光滑区域能够达到设计的7阶精度. 相比经典的7阶WENO格式, WCNS格式在精度和收敛性方面更优, 且获得相同误差时, 所耗CPU时间更短, 因此计算效率略高. 分辨率测试算例结果表明, 相比经典的7阶WENO格式, WCNS格式具有更好的分辨率和更强的间断捕捉能力. 其他一维和二维测试算例结果均表明WCNS格式精度高、分辨率高、间断捕捉能力强. 下一步工作, 进一步优化全局模板和子模板的光滑度量指标, 提高WCNS格式计算效率.
-
![]()
图 1 两种格式的计算误差与CPU时间关系曲线图
Figure 1. CPU times against numerical computing errors obtained with two schemes
![]()
图 2 基于初值条件(33)在
$t = 11$ 时刻所得数值解比较Figure 2. Comparisons of numerical solutions at time
$t = 11$ based on the initial condition (33)![]()
图 3 基于初值条件(34)在不同时刻所得数值解比较
Figure 3. Comparisons of numerical solutions at different times based on the initial condition (34)
![]()
图 4 一维(a)凸和(b)非凸Hamilton问题的数值解
Figure 4. Numerical solutions of 1D (a) convex and (b) nonconvex Hamilton problems
![]()
图 5 二维凸Hamilton问题的数值解(a)曲面图和(b)截面图
Figure 5. (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D convex Hamilton problem
![]()
图 6 二维非凸Hamilton问题的数值解(a)曲面图和(b)截面图
Figure 6. (a) Surface and (b) cross-sectional diagram of the numerical solution of the 2D nonconvex Hamilton problem
![]()
图 8 二维最优控制问题的数值解(a), (b)曲面图和(c), (d)截面图
Figure 8. (a), (b) Surfaces and (c), (d) cross-sectional diagrams of the numerical solutions of the optimal control problem
![]()
8 二维最优控制问题的数值解(a), (b)曲面图和(c), (d)截面图 (续)
8. (a), (b)Surfaces and (c), (d) cross-sectional diagrams of the numerical solutions of the optimal control problem (continued)
![]()
7 二维完全问题的数值解(a), (b) 曲面图和(c), (d)截面图
7. (a), (b) Surfaces and (c), (d) cross-sectional diagrams of numerical solutions of the fully problem
![]()
图 9 二维传播面问题的
$\varepsilon = 0$ 时数值解Figure 9. Numerical solutions of the propagating surface problem with
$\varepsilon = $ 0![]()
图 10 二维传播面问题的
$\varepsilon = 0.1$ 时数值解Figure 10. Numerical solutions of the propagating surface problem with
$\varepsilon = $ 0.1表 1 一维情形时两种格式的数值误差、精度阶和CPU运行时间
Table 1 Numerical errors, convergence rates and CPU time obtained with two schemes for the 1D case
Method $N$ ${L_1}$error Order ${L_\infty }$error Order CPU time WCNS 60 2.00 × 10−8 — 3.72 × 10−8 — 0.078125 80 2.67 × 10−9 7.00 4.90 × 10−9 7.04 0.218750 100 5.60 × 10−10 7.00 1.02 × 10−9 7.05 0.453125 120 1.56 × 10−10 7.00 2.81 × 10−10 7.05 0.796875 140 5.33 × 10−11 6.98 9.52 × 10−11 7.03 1.312500 WENO 60 7.73 × 10−8 — 5.94 × 10−7 — 0.078125 80 1.33 × 10−8 6.13 1.31 × 10−7 5.24 0.203125 100 3.36 × 10−9 6.15 4.09 × 10−8 5.23 0.453125 120 1.14 × 10−9 5.94 1.57 × 10−8 5.24 0.828125 140 4.56 × 10−10 5.93 6.97× 10−9 5.27 1.093750 
下载: 导出CSV
表 2 二维情形时两种格式的数值误差、精度阶和CPU运行时间
Table 2 Numerical errors, convergence rates and CPU time obtained with two schemes for the 2D case
Method $M \times N$ ${L_1}$error Order ${L_\infty }$error Order CPU time WCNS 60 × 60 2.95× 10−8 — 4.80× 10−8 — 1.578125 80 × 80 3.95× 10−9 6.99 6.42× 10−9 6.99 4.843750 100 × 100 8.28 × 10−10 7.00 1.34× 10−9 7.01 12.65625 120 × 120 2.31 × 10−10 7.00 3.74 × 10−10 7.02 27.85937 140 × 140 7.86 × 10−11 7.00 1.27 × 10−10 7.01 54.92187 WENO 60 × 60 2.96× 10−8 — 1.88 × 10−7 — 1.140625 80 × 80 4.48× 10−9 6.57 3.84× 10−8 5.52 3.765625 100 × 100 1.07× 10−9 6.41 1.11× 10−8 5.55 9.984375 120 × 120 3.33 × 10−10 6.41 3.94× 10−9 5.70 22.09375 140 × 140 1.24 × 10−10 6.41 1.80× 10−9 5.09 43.14062
下载: 导出CSV
-
[1] 刘瑛, 杜光勋, 全权等. 基于Hamilton-Jacobi方程的飞行器机动动作可达集分析. 自动化学报, 2016, 42(3): 347-357 (Liu Ying, Du Guangxun, Quan Quan, et al. Reachability calculation for aircraft maneuver using hamilton-jacobi function. Acta Automatica Sinica, 2016, 42(3): 347-357 (in Chinese) doi: 10.16383/j.aas.2016.c140888 [2] Chan CK, Lau KS, Zhang BL. Simulation of a premixed turbulent flame with the discrete vortex method. International Journal for Numerical Methods in Engineering, 2000, 48(4): 613-627 doi: 10.1002/(SICI)1097-0207(20000610)48:4<613::AID-NME900>3.0.CO;2-7
[3] Huang L, Wong SC, Zhang M, et al. Revisiting Hughes’ dynamic continuum model for pedestrian flow and the development of an efficient solution algorithm. Transportation Research Part B:Methodological, 2009, 43(1): 127-141 doi: 10.1016/j.trb.2008.06.003
[4] Lefèvre V, Garnica A, Lopez-Pamies O. A WENO finite-difference scheme for a new class of Hamilton-Jacobi equations in nonlinear solid mechanics. Computer Methods in Applied Mechanics and Engineering, 2019, 349(1): 17-44
[5] 陈苏婷, 李霞. 非紧空间上折现Hamilton-Jacobi方程的黏性解. 华东师范大学学报, 2022, 2022(2): 9-15 (Chen Suting, Li Xia. The viscosity solution of the discounted Hamilton-Jacobi equation in non-compact space. Proceedings of the National Academy of Sciences, 2022, 2022(2): 9-15 (in Chinese) [6] Crandall MG, Lions PL. Viscosity solutions of Hamilton-Jacobi equations. Transactions of the American Mathematical Society, 1983, 277: 1-42 doi: 10.1090/S0002-9947-1983-0690039-8
[7] Acker F, Borges RBR, Costa B. An improved WENO-Z scheme. Journal of Computational Physics, 2016, 313: 726-753 doi: 10.1016/j.jcp.2016.01.038
[8] 骆信, 吴颂平. 改进的五阶WENO-Z+ 格式. 力学学报, 2019, 51(6): 1927-1939 (Luo Xin, Wu Songping. An improved fifth-order WENO-Z + scheme. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(6): 1927-1939 (in Chinese) doi: 10.6052/0459-1879-19-249 [9] Balsara DS, Shu CW. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics, 2000, 160(2): 405-452 doi: 10.1006/jcph.2000.6443
[10] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79(1): 12-49 doi: 10.1016/0021-9991(88)90002-2
[11] Jiang GS, Peng D. Weighted ENO schemes for Hamilton-Jacobi equations. SIAM Journal on Scientific Computing, 2000, 21(6): 2126-2143 doi: 10.1137/S106482759732455X
[12] Bryson S, Levy D. High-order central WENO schemes for multidimensional Hamilton-Jacobi equations. SIAM Journal on Numerical Analysis, 2003, 41(4): 1339-1369 doi: 10.1137/S0036142902408404
[13] Qiu J, Shu CW. Hermite WENO schemes for Hamilton–Jacobi equations. Journal of Computational Physics, 2005, 204(1): 82-99 doi: 10.1016/j.jcp.2004.10.003
[14] Ha Y, Kim CH, Yang H, et al. A sixth-order weighted essentially non-oscillatory schemes based on exponential polynomials for Hamilton−Jacobi equations. Journal of Scientific Computing, 2018, 75(3): 1675-1700 doi: 10.1007/s10915-017-0603-8
[15] Cheng X, Feng J. A sixth-order finite difference WENO scheme for Hamilton−Jacobi equations. International Journal of Computer Mathematics, 2019, 96(3): 568-584 doi: 10.1080/00207160.2018.1447665
[16] Zhu J, Zheng F, Qiu J. New Finite Difference Hermite WENO Schemes for Hamilton-Jacobi Equations. Journal of Scientific Computing, 2020, 83(1): 1-21 doi: 10.1007/s10915-020-01189-x
[17] Rathan S. L1-type smoothness indicators based weighted essentially nonoscillatory scheme for Hamilton-Jacobi equations. International Journal for Numerical Methods in Fluids, 2020, 92(12): 1927-1947 doi: 10.1002/fld.4855
[18] Abedian R. A symmetrical WENO-Z scheme for solving Hamilton−Jacobi equations. International Journal of Modern Physics C, 2020, 31(3): 2050039 doi: 10.1142/S0129183120500394
[19] Kim CH, Ha Y, Yang H, et al. A third-order WENO scheme based on exponential polynomials for Hamilton-Jacobi equations. Applied Numerical Mathematics, 2021, 165: 167-183 doi: 10.1016/j.apnum.2021.01.020
[20] Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes. Journal of Computational Physics, 2000, 165(1): 22-44 doi: 10.1006/jcph.2000.6594
[21] Deng X, Maekawa H. Compact high-order accurate nonlinear schemes. Journal of Computational Physics, 1997, 130(1): 77-91 doi: 10.1006/jcph.1996.5553
[22] Nonomura T, Fujii K. Effects of difference scheme type in high-order weighted compact nonlinear schemes. Journal of Computational Physics, 2009, 228(10): 3533-3539 doi: 10.1016/j.jcp.2009.02.018
[23] Zhang S, Jiang S, Shu CW. Development of nonlinear weighted compact schemes with increasingly higher order accuracy. Journal of Computational Physics, 2008, 227(15): 7294-7321 doi: 10.1016/j.jcp.2008.04.012
[24] Deng X, Jiang Y, Mao M, et al. A family of hybrid cell-edge and cell-node dissipative compact schemes satisfying geometric conservation law. Computers & Fluids, 2015, 116: 29-45
[25] Nonomura T, Fujii K. Robust explicit formulation of weighted compact nonlinear scheme. Computers & Fluids, 2013, 85: 8-18
[26] Wong ML, Lele SK. High-order localized dissipation weighted compact nonlinear scheme for shock-and interface-capturing in compressible flows. Journal of Computational Physics, 2017, 339: 179-209 doi: 10.1016/j.jcp.2017.03.008
[27] Zhao G, Sun M, Xie S, et al. Numerical dissipation control in an adaptive WCNS with a new smoothness indicator. Applied Mathematics and Computation, 2018, 330: 239-253 doi: 10.1016/j.amc.2018.01.019
[28] 洪正, 叶正寅. 各向同性湍流通过正激波的演化特征研究. 力学学报, 2018, 50(6): 1356-1367 (Hong Zheng, Ye Zhengyin. Study on evolution characteristics of isotropic turbulence passing through a normal shock wave. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(6): 1356-1367 (in Chinese) doi: 10.6052/0459-1879-18-129 [29] Jiang YQ, Zhou SG, Zhang X, et al. High-order weighted compact nonlinear scheme for one-and two-dimensional Hamilton-Jacobi equations. Applied Numerical Mathematics, 2022, 171: 353-368 doi: 10.1016/j.apnum.2021.09.012
[30] Hiejima T. A high-order weighted compact nonlinear scheme for compressible flows. Computers & Fluids, 2022, 232: 105199
[31] Jiang YQ, Zhou SG, Zhang X, et al. High order all-speed semi-implicit weighted compact nonlinear scheme for the isentropic Navier–Stokes equations. Journal of Computational and Applied Mathematics, 2022, 411: 114272 doi: 10.1016/j.cam.2022.114272
[32] Ma Y, Yan ZG, Liu H, et al. Improved weighted compact nonlinear scheme for implicit large-eddy simulations. Computers & Fluids, 2022, 240: 105412
[33] Gottlieb S, Shu CW, Tadmor E. Strong stability-preserving high-order time discretization methods. SIAM Review, 2001, 43(1): 89-112 doi: 10.1137/S003614450036757X
-
期刊类型引用(6)
1. 都凯，任彦强，侯勇，李小强，陈帅峰，孙亮，左晓姣，袁晓光. 精确预测近平面应变载荷下各向异性行为的屈服准则参数识别策略. 中国机械工程. 2023(19): 2381-2393 . 
百度学术
2. 王海波，肖旭，阎昱，李强，何东. 流动法则和屈服准则对板料成形极限预测的影响. 塑性工程学报. 2022(03): 157-165 . 百度学术
3. 方刚，陈祝，雷丽萍. 非关联本构模型在铝合金板料成形有限元模拟中的应用. 塑性工程学报. 2021(06): 8-18 . 百度学术
4. 贾然，赵桂平. 泡沫铝本构行为研究进展. 力学学报. 2020(03): 603-622 . 本站查看
5. 陶文斌，陈铁林，张群，王泽军，王荣鑫. 考虑非均匀应力场巷道弹塑性分区的锚杆力学分析. 吉林大学学报(工学版). 2020(06): 2131-2140 . 百度学术
6. Haibo Wang，Mingliang Men，Yu Yan，Min Wan，Qiang Li. Prediction of Eight Earings in Deep Drawing of 5754O Aluminum Alloy Sheet. Chinese Journal of Mechanical Engineering. 2019(05): 146-154 . 必应学术
其他类型引用(5)
计量
- 文章访问数: 764
- HTML全文浏览量: 212
- PDF下载量: 108
- 被引次数: 11
