高分辨率格式理论
一个核心概念：人工粘性
考虑经典的双曲守恒律方程
∂u∂t+∂f∂x=0{{\partial u} \over {\partial t}} + {{\partial f} \over {\partial x}} = 0∂t∂u+∂x∂f=0
可以写成守恒形式的数值格式
uin+1=uin−λ(f^i+1/2n−f^i+1/2n)u_i^{n + 1} = u_i^n - \lambda \left( {\hat f_{i + 1/2}^n - \hat f_{i + 1/2}^n} \right)uin+1=uin−λ(f^i+1/2n−f^i+1/2n)
式中λ\lambdaλ是时间步长与空间步长之比。f^{\hat f}f^则是面上的数值通量。各种格式的最终目的就是设法给出数值通量的近似。
为了减小间断处的震荡，一种有效的方法是在方程右端添加人工粘性项
∂u∂t+∂f∂x=Δx2Δt∂∂x(Q∂u∂x){{\partial u} \over {\partial t}} + {{\partial f} \over {\partial x}} = {{\Delta {x^2}} \over {\Delta t}}{\partial \over {\partial x}}\left( {Q{{\partial u} \over {\partial x}}} \right)∂t∂u+∂x∂f=ΔtΔx2∂x∂(Q∂x∂u)相似地，也可差分得到下式
uin+1=uin−λ(f^i+1/2n−f^i+1/2n)+1/2(Qi+1/2Δ+uin−Qi−1/2Δ−uin)u_i^{n + 1} = u_i^n - \lambda \left( {\hat f_{i + 1/2}^n - \hat f_{i + 1/2}^n} \right) +1/2 \left( {{Q_{i + 1/2}}{\Delta ^ + }u_i^n - {Q_{i - 1/2}}{\Delta ^ - }u_i^n} \right)uin+1=uin−λ(f^i+1/2n−f^i+1/2n)+1/2(Qi+1/2Δ+uin−Qi−1/2Δ−uin)
可以将右侧第三项与数值通量合并，得到新的数值通量
fˉi+1/2n=f^i+1/2n−Qi+1/22λΔ+uin\bar f_{i + 1/2}^n = \hat f_{i + 1/2}^n - {{{Q_{i + 1/2}}} \over {2\lambda }}{\Delta ^ + }u_i^nfˉi+1/2n=f^i+1/2n−2λQi+1/2Δ+uinfˉi−1/2n=f^i−1/2n−Qi−1/22λΔ−uin\bar f_{i - 1/2}^n = \hat f_{i - 1/2}^n - {{{Q_{i - 1/2}}} \over 2\lambda }{\Delta ^ - }u_i^nfˉi−1/2n=f^i−1/2n−2λQi−1/2Δ−uin
这里可以看出“人工粘性”可增加逆梯度方向上的通量，即如果ui+1n>uinu_{i + 1}^n > u_i^nui+1n>uin，那么“人工粘性”将会减少i+1/2i + 1/2i+1/2面上的正向通量。相反地，负人工粘性则会增加i+1/2i + 1/2i+1/2面上的通量。合适地选取、改变、调整负人工粘性则是开发高分辨率格式的核心问题。
经典格式
中心差分格式
f^i+1/2n=fi+fi+12\hat f_{i + 1/2}^n = {{{f_i} + {f_{i + 1}}} \over 2}f^i+1/2n=2fi+fi+1
一阶迎风格式
f^i+1/2n=fi+fi+12−12∣ai+1/2n∣(ui+1−ui)\hat f_{i + 1/2}^n = {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\left| {a_{i + 1/2}^n} \right|\left( {{u_{i + 1}} - {u_i}} \right)f^i+1/2n=2fi+fi+1−21∣∣∣ai+1/2n∣∣∣(ui+1−ui)
式中右端第二项为一阶迎风格式引入的数值粘性
Lax-Fridrichs一阶格式
f^i+1/2n=fi+fi+12−12λ(ui+1−ui)\hat f_{i + 1/2}^n = {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over {2\lambda }}\left( {{u_{i + 1}} - {u_i}} \right)f^i+1/2n=2fi+fi+1−2λ1(ui+1−ui)
第二项是数值粘性，该格式的数值粘性非常大
Lax-Wendroff二阶格式
f^i+1/2n=fi+fi+12−12λ(ai+1/2n)2(ui+1−ui)\hat f_{i + 1/2}^n = {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\lambda {\left( {a_{i + 1/2}^n} \right)^2}\left( {{u_{i + 1}} - {u_i}} \right)f^i+1/2n=2fi+fi+1−21λ(ai+1/2n)2(ui+1−ui)
注意："数值粘性"来源于格式本身，"人工粘性"来源于人为添加的扩散项，虽然表象不同，但二者的本质是相同的。
高分辨率TVD格式的原理
尝试在1阶TVD格式的基础上添加负的"人工粘性“。这个“人工粘性”必须是非线性的，使得：在光滑区域达到二阶格式，而在间断区域回归到一阶格式。
A 通量修正法
Harten1提出的通量修正法即是在1阶TVD格式的基础上修正物理通量fff以部分抵消原格式的截断误差。这个修正量（事实上的“人工粘性”）必须是非线性的，使得：在光滑区域完全抵消低阶格式的数值粘性，并趋近于二阶L-W格式的数值粘性，达到二阶格式。而在间断区域回归到一阶格式。为了实现这一目标，必须平衡好fˉi+1/2n\bar f_{i + 1/2}^nfˉi+1/2n中蕴含的“数值粘性”部分与“人工粘性”部分。
B 通量限制器
与Harten的想法类似，但是不再修正"物理通量fff"，而是直接修正"数值通量fˉi+1/2n\bar f_{i + 1/2}^nfˉi+1/2n"。将高阶格式的数值通量写为低阶格式的数值通量与反扩散通量之和，为了使格式总变差不增（TVD条件），需要限制反扩散通量的大小。这就是通量限制器。反扩散通量实质上就是负的"人工粘性"。
比如，取使用L-W格式作为高阶格式
[f^i+1/2n]High=fi+fi+12−12λ(ai+1/2n)2(ui+1−ui){\left[ {\hat f_{i + 1/2}^n} \right]_{{\rm{High}}}} = {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\lambda {\left( {a_{i + 1/2}^n} \right)^2}\left( {{u_{i + 1}} - {u_i}} \right)[f^i+1/2n]High=2fi+fi+1−21λ(ai+1/2n)2(ui+1−ui)
再使用一阶迎风格式作为低价格式
[f^i+1/2n]Low=fi+fi+12−12∣ai+1/2n∣(ui+1−ui){\left[ {\hat f_{i + 1/2}^n} \right]_{{\rm{Low}}}} = {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\left| {a_{i + 1/2}^n} \right|\left( {{u_{i + 1}} - {u_i}} \right)[f^i+1/2n]Low=2fi+fi+1−21∣∣∣ai+1/2n∣∣∣(ui+1−ui)
两种格式进行组合
f^i+1/2n=[f^i+1/2n]Low+φ([f^i+1/2n]High−[f^i+1/2n]Low)\hat f_{i + 1/2}^n = {\left[ {\hat f_{i + 1/2}^n} \right]_{{\rm{Low}}}} + \varphi \left( {{{\left[ {\hat f_{i + 1/2}^n} \right]}_{{\rm{High}}}} - {{\left[ {\hat f_{i + 1/2}^n} \right]}_{{\rm{Low}}}}} \right)f^i+1/2n=[f^i+1/2n]Low+φ([f^i+1/2n]High−[f^i+1/2n]Low)=fi+fi+12−12∣ai+1/2n∣(ui+1−ui)+φ(−12λ(ai+1/2n)2(ui+1−ui)+12∣ai+1/2n∣(ui+1−ui))= {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\left| {a_{i + 1/2}^n} \right|\left( {{u_{i + 1}} - {u_i}} \right) + \varphi \left( { - {1 \over 2}\lambda {{\left( {a_{i + 1/2}^n} \right)}^2}\left( {{u_{i + 1}} - {u_i}} \right) + {1 \over 2}\left| {a_{i + 1/2}^n} \right|\left( {{u_{i + 1}} - {u_i}} \right)} \right)=2fi+fi+1−21∣∣∣ai+1/2n∣∣∣(ui+1−ui)+φ(−21λ(ai+1/2n)2(ui+1−ui)+21∣∣∣ai+1/2n∣∣∣(ui+1−ui))=fi+fi+12−12∣ai+1/2n∣(ui+1−ui)+φsgn(νi+1/2n)−νi+1/2n2ai+1/2n(ui+1−ui)= {{{f_i} + {f_{i + 1}}} \over 2} - {1 \over 2}\left| {a_{i + 1/2}^n} \right|\left( {{u_{i + 1}} - {u_i}} \right) + \varphi {{{\mathop{\rm sgn}} \left( {\nu _{i + 1/2}^n} \right) - \nu _{i + 1/2}^n} \over 2}a_{i + 1/2}^n\left( {{u_{i + 1}} - {u_i}} \right)=2fi+fi+1−21∣∣∣ai+1/2n∣∣∣(ui+1−ui)+φ2sgn(νi+1/2n)−νi+1/2nai+1/2n(ui+1−ui)
式中的φ\varphiφ表示通量限制器，对于光滑区域取1，达到二阶L-W格式，对于间断区域取0，恢复到一阶迎风格式。首先使用变量rrr描述解的光滑程度
ri+1/2={ui−ui−1ui+1−ui,ai+1/2>0ui+2−ui+1ui+1−ui,ai+1/2<0r_{i+1 / 2}=\left\{\begin{array}{l} \frac{u_{i}-u_{i-1}}{u_{i+1}-u_{i}}, a_{i+1 / 2}>0 \\ \frac{u_{i+2}-u_{i+1}}{u_{i+1}-u_{i}}, a_{i+1 / 2}<0 \end{array}\right. ri+1/2={ui+1−uiui−ui−1,ai+1/2>0ui+1−uiui+2−ui+1,ai+1/2<0
然后使用rrr描述通量限制器φ\varphiφ，通量限制器的表达式有许多，详见Sweby2的工作
C 斜率限制器
采用分片线性函数重构解，在单元面上应用TVD模板。
f^i+1/2(ui,ui+1)→f^i+1/2(ui+1/2L,ui+1/2R){\hat f_{i + 1/2}}\left( {{u_i},{u_{i + 1}}} \right) \to {\hat f_{i + 1/2}}\left( {u_{i + 1/2}^L,u_{i + 1/2}^R} \right)f^i+1/2(ui,ui+1)→f^i+1/2(ui+1/2L,ui+1/2R)ui+1/2L=ui+siΔx2,ui+1/2R=ui+1−si+1Δx2u_{i + 1/2}^L = {u_i} + {s_i}{{\Delta x} \over 2},u_{i + 1/2}^R = {u_{i + 1}} - {s_{i + 1}}{{\Delta x} \over 2}ui+1/2L=ui+si2Δx,ui+1/2R=ui+1−si+12Δx可以证明，在一定条件下，前者是一阶TVD的，那么后者是二阶TVD的
需要使用限制器来限制线性函数的斜率以实现TVD格式。其本质与通量限制器相同。限制后的斜率与通量限制器之间的关系是
si=(ui+1−uiΔx)φ(ri+1/2){s _{i}} = \left( {{{{u_{i + 1}} - {u_i}} \over {\Delta x}}} \right)\varphi \left( {{r_{i + 1/2}}} \right)si=(Δxui+1−ui)φ(ri+1/2)具体可参考MUSCL格式
HARTEN A. High resolution schemes for hyperbolic conservation laws[J]. Journal of Computationalphysics, 1983(49): 357-393. ↩︎
SWEBY P K. High resolution schemes using flux limiters for hyperbolic conservation laws[J]. Siam Journal On Numerical Analysis, 1984, 21(5): 995-1011. ↩︎