Kurganov-Tadmor二阶中心格式：理论介绍

简介

CFD的核心问题是求解双曲偏微分方程
$$\frac{\partial }{\partial t}u(x,t)+\frac{\partial }{\partial x}f(u(x,t))=0$$在CFD中，双曲偏微分方程一般使用Godunov型迎风格式求解。但是这种迎风格式往往实施起来比较复杂（需要特征分解），如果能使用中心格式实现离散，则可以简化编程、提高计算效率。

Kurganov-Tadmor格式¹（简称KT）可以实现二阶精度，是一种中心离散格式。根据此格式openFoam开发了rhoCentralFoam求解器，以求解可压缩流。本文仅介绍格式的推导思路，具体在openFoam中的实现方法将由后续文章给出。

KT格式有两种形式，它们的推导方式如下：
（1）全离散形式由分片线性解积分所得。
（2）半离散形式由Rusanov离散格式演化而来。

全离散格式

使用分片线性函数完成解的重构，重构后的解为
$$\tilde{u}\left(x,t^n\right):=\sum _j\left[{\bar{u}}_j^n+{\left(u_x\right)}_j^n\left(x-x_j\right)\right]1_{\left[x_{j-1/2},x_{j+1/2}\right]},x_{j\pm 1/2}:=x_j\pm \frac{\Delta x}{2}$$
图中$x_{j+\frac{1}{2},l}^n$表示在一个时间步后，扰动由$x_{j+\frac{1}{2}}$向左传播到的位置。所以有
$$x_{j+\frac{1}{2}}-x_{j+\frac{1}{2},l}^n=a_{j+\frac{1}{2}}^n\Delta t$$同理有
$$x_{j+\frac{1}{2},r}^n-x_{j+\frac{1}{2}}=a_{j+\frac{1}{2}}^n\Delta t$$式中$a_{j+\frac{1}{2}}^n$是界面处扰动的传播速度。可使用朗金-雨果纽关系式确定。

(1) 在$\left[x_{j+\frac{1}{2},l}^n,x_{j+\frac{1}{2},r}^n\right]\times \left[t^n,t^{n+1}\right]$范围内，对双曲守恒方程积分可得

$$\begin{array}{rl}{w}_{j+1/2}^{n+1}=& \frac{u_j^n+u_{j+1}^n}{2}+\frac{\Delta x-a_{j+1/2}^n\Delta t}{4}\left({\left(u_x\right)}_j^n-{\left(u_x\right)}_{j+1}^n\right)\\ & -\frac{1}{2a_{j+1/2}^n}\left[f\left(u_{j+1/2,r}^{n+1/2}\right)-f\left(u_{j+1/2,l}^{n+1/2}\right)\right],\end{array}$$
式中，$w_{j+1/2}^{n+1}$表示$\left[x_{j+\frac{1}{2},l}^n,x_{j+\frac{1}{2},r}^n\right]$范围内的平均积分结果，上角标$n+1/2$表示半时间步，可用泰勒展开近似
$$u_{j+1/2,l}^{n+1/2}:=u_{j+1/2,l}^n-\frac{\Delta t}{2}f{\left(u_{j+1/2,l}^n\right)}_x,u_{j+1/2,l}^n:=u_j^n+\Delta x{\left(u_x\right)}_j^n\left(\frac{1}{2}-\lambda a_{j+1/2}^n\right)$$同理有
$$u_{j+1/2,r}^{n+1/2}:=u_{j+1/2,r}^n-\frac{\Delta t}{2}f{\left(u_{j+1/2,r}^n\right)}_x,u_{j+1/2,r}^n:=u_{j+1}^n-\Delta x{\left(u_x\right)}_{j+1}^n\left(\frac{1}{2}-\lambda a_{j+1/2}^n\right)$$

(2) 在$\left[x_{j-\frac{1}{2},r}^n,x_{j+\frac{1}{2},l}^n\right]\times \left[t^n,t^{n+1}\right]$范围内，对双曲守恒方程积分可得
$$\begin{array}{rl}{w}_j^{n+1}=& u_j^n+\frac{\Delta t}{2}\left(a_{j-1/2}^n-a_{j+1/2}^n\right){\left(u_x\right)}_j^n\\ & -\frac{\lambda }{1-\lambda \left(a_{j-1/2}^n+a_{j+1/2}^n\right)}\left[f\left(u_{j+1/2,l}^{n+1/2}\right)-f\left(u_{j-1/2,r}^{n+1/2}\right)\right]\end{array}$$

综上得到了区间$\left[x_{j+\frac{1}{2},l}^n,x_{j+\frac{1}{2},r}^n\right]$与区间$\left[x_{j-\frac{1}{2},r}^n,x_{j+\frac{1}{2},l}^n\right]$上解的均值。

为了得到最终的$u_j^{n+1}$,需要将区间$\left[x_{j-\frac{1}{2}}^n,x_{j-\frac{1}{2},r}^n\right]$、区间$\left[x_{j-\frac{1}{2},r}^n,x_{j+\frac{1}{2},l}^n\right]$和$\left[x_{j+\frac{1}{2},l}^n,x_{j+\frac{1}{2}}^n\right]$的解积分平均，最终得到
$$\begin{array}{rl}{u}_j^{n+1}=& \frac{1}{\Delta x}{\int }_{x_{j-1/2}}^{x_{j+1/2}}\tilde{w}\left(\xi ,t^{n+1}\right)d\xi =\lambda a_{j-1/2}^n{w}_{j-1/2}^{n+1}+\left[1-\lambda \left(a_{j-1/2}^n+a_{j+1/2}^n\right)\right]w_j^{n+1}\\ & +\lambda a_{j+1/2}^n{w}_{j+1/2}^{n+1}+\frac{\Delta x}{2}\left[{\left(\lambda a_{j-1/2}^n\right)}^2{\left(u_x\right)}_{j-1/2}^{n+1}-{\left(\lambda a_{j+1/2}^n\right)}^2{\left(u_x\right)}_{j+1/2}^{n+1}\right]\end{array}$$
这就是全离散K-T格式，具体的minmod斜率限制器可参见论文原文。斜率限制器的目的是使得此二阶格式TVD，在间断解处退回一阶格式。

半离散格式

半离散格式的推导相对简单，首先考虑一阶Rusanov格式
$$u_j^{n+1}=u_j^n-\frac{\lambda }{2}\left[f\left(u_{j+1}^n\right)-f\left(u_{j-1}^n\right)\right]+\frac{1}{2}\left[\lambda a_{j+1/2}^n\left(u_{j+1}^n-u_j^n\right)-\lambda a_{j-1/2}^n\left(u_j^n-u_{j-1}^n\right)\right]$$
将其写为半离散格式（左端项是时间导数）
$$\frac{\partial u}{\partial t}=-\frac{1}{\Delta x}\left(\frac{f\left(u_{j+1}^n\right)-f\left(u_{j-1}^n\right)}{2}-\frac{a_{j+1/2}^n\left(u_{j+1}^n-u_j^n\right)-a_{j-1/2}^n\left(u_j^n-u_{j-1}^n\right)}{2}\right)$$
进一步写为守恒形式
$$\frac{\partial u}{\partial t}=-\frac{1}{\Delta x}\left(H_{j+\frac{1}{2}}-H_{j-\frac{1}{2}}\right)$$其中
$$H_{j+\frac{1}{2}}=\frac{f\left(u_{j+1}^n\right)+f\left(u_j^n\right)}{2}-\frac{a_{j+1/2}^n}{2}\left(u_{j+1}^n-u_j^n\right)$$$$H_{j-\frac{1}{2}}=\frac{f\left(u_j^n\right)+f\left(u_{j-1}^n\right)}{2}-\frac{a_{j-1/2}^n}{2}\left(u_j^n-u_{j-1}^n\right)$$
为了将格式提升为二阶，则使用分段线性将解重构，有
$$H_{j+1/2}(t):=\frac{f\left(u_{j+1/2}^{+}(t)\right)+f\left(u_{j+1/2}^{-}(t)\right)}{2}-\frac{a_{j+1/2}(t)}{2}\left[u_{j+1/2}^{+}(t)-u_{j+1/2}^{-}(t)\right]$$
式中
$$u_{j+1/2}^{+}:=u_{j+1}(t)-\frac{\Delta x}{2}{\left(u_x\right)}_{j+1}(t),u_{j+1/2}^{-}:=u_j(t)+\frac{\Delta x}{2}{\left(u_x\right)}_j(t)$$
$H_{j-1/2}$的计算公式不再赘述。

这就是半离散K-T格式，具体的minmod斜率限制器可参见论文原文。斜率限制器的目的是使得此二阶格式TVD，在间断解处退回一阶格式。

半离散格式只完成了空间离散，时间离散可使用TVD-龙格库塔法实现。

拓展阅读

为了进一步学习Kurganov-Tadmor以及其应用，有如下几篇关键文献：
（1）2001年，文章SEMIDISCRETE CENTRAL-UPWIND SCHEMES FOR HYPERBOLIC CONSERVATION LAWS AND HAMILTON–JACOBI EQUATIONS拓展了KT格式，提出了中心迎风格式。该格式也是openFoam所采用的格式。
（2）文章Central-upwind schemes for hyperbolic conservation laws简略地总结了中心迎风格式。
（3）2010年，文章Implementation of semi-discrete, non-staggered central schemes in a colocated, polyhedral, finite volume framework, for highspeed viscous flows将中心迎风格式引入openFoam,开发了rhoCentralFoam求解器
（4）2016年，文章Density-based solver for all Mach number flows提出了适用于全马赫数的openFoam密度基求解器
（5）2017年，文章A hybrid pressure-based solver for nonideal single-phase fluid flows at all speeds提出了适用于全马赫数的openFoam压力基求解器

---

1. KURGANOV A, TADMOR E. New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations[J]. Journal of Computational Physics, 2000, 160(1): 241-282. ↩︎