Kurganov-Tadmor二阶中心格式：理论介绍
简介
CFD的核心问题是求解双曲偏微分方程
∂∂tu(x,t)+∂∂xf(u(x,t))=0\frac{\partial}{\partial t} u(x, t)+\frac{\partial}{\partial x} f(u(x, t))=0 ∂t∂u(x,t)+∂x∂f(u(x,t))=0在CFD中，双曲偏微分方程一般使用Godunov型迎风格式求解。但是这种迎风格式往往实施起来比较复杂（需要特征分解），如果能使用中心格式实现离散，则可以简化编程、提高计算效率。
Kurganov-Tadmor格式1（简称KT）可以实现二阶精度，是一种中心离散格式。根据此格式openFoam开发了rhoCentralFoam求解器，以求解可压缩流。本文仅介绍格式的推导思路，具体在openFoam中的实现方法将由后续文章给出。
KT格式有两种形式，它们的推导方式如下：
（1）全离散形式由分片线性解积分所得。
（2）半离散形式由Rusanov离散格式演化而来。
全离散格式
使用分片线性函数完成解的重构，重构后的解为
u~(x,tn):=∑j[uˉjn+(ux)jn(x−xj)]1[xj−1/2,xj+1/2],xj±1/2:=xj±Δx2\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] \mathbf{1}_{\left[x_{j-1 / 2}, x_{j+1 / 2}\right]}, \quad x_{j \pm 1 / 2}:=x_{j} \pm \frac{\Delta x}{2} u~(x,tn):=j∑[uˉjn+(ux)jn(x−xj)]1[xj−1/2,xj+1/2],xj±1/2:=xj±2Δx
图中xj+12,lnx_{j + {1 \over 2},l}^nxj+21,ln表示在一个时间步后，扰动由xj+12{x_{j + {1 \over 2}}}xj+21向左传播到的位置。所以有
xj+12−xj+12,ln=aj+12nΔt{x_{j + {1 \over 2}}} - x_{j + {1 \over 2},l}^n = a_{j + {1 \over 2}}^n\Delta txj+21−xj+21,ln=aj+21nΔt同理有
xj+12,rn−xj+12=aj+12nΔtx_{j + {1 \over 2},r}^n - {x_{j + {1 \over 2}}} = a_{j + {1 \over 2}}^n\Delta txj+21,rn−xj+21=aj+21nΔt式中aj+12na_{j + {1 \over 2}}^naj+21n是界面处扰动的传播速度。可使用朗金-雨果纽关系式确定。
(1) 在[xj+12,ln,xj+12,rn]×[tn,tn+1]\left[ {x_{j + {1 \over 2},l}^n,x_{j + {1 \over 2},r}^n} \right] \times \left[ {{t^n},{t^{n + 1}}} \right][xj+21,ln,xj+21,rn]×[tn,tn+1]范围内，对双曲守恒方程积分可得
wj+1/2n+1=ujn+uj+1n2+Δx−aj+1/2nΔt4((ux)jn−(ux)j+1n)−12aj+1/2n[f(uj+1/2,rn+1/2)−f(uj+1/2,ln+1/2)],\begin{aligned} 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}{2 a_{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{aligned} wj+1/2n+1=2ujn+uj+1n+4Δx−aj+1/2nΔt((ux)jn−(ux)j+1n)−2aj+1/2n1[f(uj+1/2,rn+1/2)−f(uj+1/2,ln+1/2)],
式中，wj+1/2n+1w_{j+1 / 2}^{n+1}wj+1/2n+1表示[xj+12,ln,xj+12,rn]\left[ {x_{j + {1 \over 2},l}^n,x_{j + {1 \over 2},r}^n} \right][xj+21,ln,xj+21,rn]范围内的平均积分结果，上角标n+1/2n+1/2n+1/2表示半时间步，可用泰勒展开近似
uj+1/2,ln+1/2:=uj+1/2,ln−Δt2f(uj+1/2,ln)x,uj+1/2,ln:=ujn+Δx(ux)jn(12−λaj+1/2n)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}, \quad 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)uj+1/2,ln+1/2:=uj+1/2,ln−2Δtf(uj+1/2,ln)x,uj+1/2,ln:=ujn+Δx(ux)jn(21−λaj+1/2n)同理有
uj+1/2,rn+1/2:=uj+1/2,rn−Δt2f(uj+1/2,rn)x,uj+1/2,rn:=uj+1n−Δx(ux)j+1n(12−λaj+1/2n)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}, \quad 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) uj+1/2,rn+1/2:=uj+1/2,rn−2Δtf(uj+1/2,rn)x,uj+1/2,rn:=uj+1n−Δx(ux)j+1n(21−λaj+1/2n)
(2) 在[xj−12,rn,xj+12,ln]×[tn,tn+1]\left[ {x_{j - {1 \over 2},r}^n,x_{j + {1 \over 2},l}^n} \right] \times \left[ {{t^n},{t^{n + 1}}} \right][xj−21,rn,xj+21,ln]×[tn,tn+1]范围内，对双曲守恒方程积分可得
wjn+1=ujn+Δt2(aj−1/2n−aj+1/2n)(ux)jn−λ1−λ(aj−1/2n+aj+1/2n)[f(uj+1/2,ln+1/2)−f(uj−1/2,rn+1/2)]\begin{aligned} 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{aligned} wjn+1=ujn+2Δt(aj−1/2n−aj+1/2n)(ux)jn−1−λ(aj−1/2n+aj+1/2n)λ[f(uj+1/2,ln+1/2)−f(uj−1/2,rn+1/2)]
综上得到了区间[xj+12,ln,xj+12,rn]\left[ {x_{j + {1 \over 2},l}^n,x_{j + {1 \over 2},r}^n} \right][xj+21,ln,xj+21,rn]与区间[xj−12,rn,xj+12,ln]\left[ {x_{j - {1 \over 2},r}^n,x_{j + {1 \over 2},l}^n} \right][xj−21,rn,xj+21,ln]上解的均值。
为了得到最终的ujn+1u_j^{n + 1}ujn+1,需要将区间[xj−12n,xj−12,rn]\left[ {x_{j - {1 \over 2}}^n,x_{j - {1 \over 2},r}^n} \right][xj−21n,xj−21,rn]、区间[xj−12,rn,xj+12,ln]\left[ {x_{j - {1 \over 2},r}^n,x_{j + {1 \over 2},l}^n} \right][xj−21,rn,xj+21,ln]和[xj+12,ln,xj+12n]\left[ {x_{j + {1 \over 2},l}^n,x_{j + {1 \over 2}}^n} \right][xj+21,ln,xj+21n]的解积分平均，最终得到
ujn+1=1Δx∫xj−1/2xj+1/2w~(ξ,tn+1)dξ=λaj−1/2nwj−1/2n+1+[1−λ(aj−1/2n+aj+1/2n)]wjn+1+λaj+1/2nwj+1/2n+1+Δx2[(λaj−1/2n)2(ux)j−1/2n+1−(λaj+1/2n)2(ux)j+1/2n+1]\begin{aligned} 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{aligned} ujn+1=Δx1∫xj−1/2xj+1/2w~(ξ,tn+1)dξ=λaj−1/2nwj−1/2n+1+[1−λ(aj−1/2n+aj+1/2n)]wjn+1+λaj+1/2nwj+1/2n+1+2Δx[(λaj−1/2n)2(ux)j−1/2n+1−(λaj+1/2n)2(ux)j+1/2n+1]
这就是全离散K-T格式，具体的minmod斜率限制器可参见论文原文。斜率限制器的目的是使得此二阶格式TVD，在间断解处退回一阶格式。
半离散格式
半离散格式的推导相对简单，首先考虑一阶Rusanov格式
ujn+1=ujn−λ2[f(uj+1n)−f(uj−1n)]+12[λaj+1/2n(uj+1n−ujn)−λaj−1/2n(ujn−uj−1n)]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] ujn+1=ujn−2λ[f(uj+1n)−f(uj−1n)]+21[λaj+1/2n(uj+1n−ujn)−λaj−1/2n(ujn−uj−1n)]
将其写为半离散格式（左端项是时间导数）
∂u∂t=−1Δx(f(uj+1n)−f(uj−1n)2−aj+1/2n(uj+1n−ujn)−aj−1/2n(ujn−uj−1n)2){{\partial u} \over {\partial t}} = - {1 \over {\Delta x}}\left( {{{f\left( {u_{j + 1}^n} \right) - f\left( {u_{j - 1}^n} \right)} \over 2} - {{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)} \over 2}} \right)∂t∂u=−Δx1(2f(uj+1n)−f(uj−1n)−2aj+1/2n(uj+1n−ujn)−aj−1/2n(ujn−uj−1n))
进一步写为守恒形式
∂u∂t=−1Δx(Hj+12−Hj−12){{\partial u} \over {\partial t}} = - {1 \over {\Delta x}}\left( {{H_{j + {1 \over 2}}} - {H_{j - {1 \over 2}}}} \right)∂t∂u=−Δx1(Hj+21−Hj−21)其中
Hj+12=f(uj+1n)+f(ujn)2−aj+1/2n2(uj+1n−ujn){H_{j + {1 \over 2}}} = {{f\left( {u_{j + 1}^n} \right) + f\left( {u_j^n} \right)} \over 2} - {{a_{j + 1/2}^n} \over 2}\left( {u_{j + 1}^n - u_j^n} \right)Hj+21=2f(uj+1n)+f(ujn)−2aj+1/2n(uj+1n−ujn)Hj−12=f(ujn)+f(uj−1n)2−aj−1/2n2(ujn−uj−1n){H_{j - {1 \over 2}}} = {{f\left( {u_j^n} \right) + f\left( {u_{j - 1}^n} \right)} \over 2} - {{a_{j - 1/2}^n} \over 2}\left( {u_j^n - u_{j - 1}^n} \right)Hj−21=2f(ujn)+f(uj−1n)−2aj−1/2n(ujn−uj−1n)
为了将格式提升为二阶，则使用分段线性将解重构，有
Hj+1/2(t):=f(uj+1/2+(t))+f(uj+1/2−(t))2−aj+1/2(t)2[uj+1/2+(t)−uj+1/2−(t)]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] Hj+1/2(t):=2f(uj+1/2+(t))+f(uj+1/2−(t))−2aj+1/2(t)[uj+1/2+(t)−uj+1/2−(t)]
式中
uj+1/2+:=uj+1(t)−Δx2(ux)j+1(t),uj+1/2−:=uj(t)+Δx2(ux)j(t)u_{j+1 / 2}^{+}:=u_{j+1}(t)-\frac{\Delta x}{2}\left(u_{x}\right)_{j+1}(t), \quad u_{j+1 / 2}^{-}:=u_{j}(t)+\frac{\Delta x}{2}\left(u_{x}\right)_{j}(t) uj+1/2+:=uj+1(t)−2Δx(ux)j+1(t),uj+1/2−:=uj(t)+2Δx(ux)j(t)
Hj−1/2H_{j-1 / 2}Hj−1/2的计算公式不再赘述。
这就是半离散K-T格式，具体的minmod斜率限制器可参见论文原文。斜率限制器的目的是使得此二阶格式TVD，在间断解处退回一阶格式。
半离散格式只完成了空间离散，时间离散可使用TVD-龙格库塔法实现。
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. ↩︎