FVM in CFD 学习笔记_第12章_高分辨率格式

2020年06月10日 阅读数:241
这篇文章主要向大家介绍FVM in CFD 学习笔记_第12章_高分辨率格式,主要内容包括基础应用、实用技巧、原理机制等方面,希望对大家有所帮助。

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 12 High Resolution Schemeshtml

本章继续讲解对流项格式的发展,讨论如何对高阶(High Order HO)格式施加有界性来产生高分辨率(High Resolution HR)格式。HR格式是HO廓线和对流有界准则(Convection Boundness Criterion CBC)的结合,以保证解不会出现振荡。将介绍发展HR格式的归一化变量公式(Normalized Variable Formulation(NVF))和总变差衰减(Total Variation Diminishing(TVD))框架,这俩方法看起来不一样但实际上几乎是同样的。分别展现了在NVF和TVD中用于使HR格式更具形象化的归一化变量图(Normalized Variable Diagram (NVD))和Sweby(或 r ψ r-\psi )图。针对NVF和TVD,具体给出了一些HR格式的函数关系。还有上一章讲到的延迟修正(Deferred Correction(DC)),介绍了两种技术来实如今结构网格和非结构网格的HO和HR格式,即,背风加权因子(Downwind Weighing Factor(DWF))方法和归一化加权因子(Normalized Weighing Factor(NWF))方法。web

1 归一化变量公式(Normalized Variable Formulation(NVF))

归一化变量公式(NVF)是描述和分析高分辨率(HR)格式的一个框架,其由Leonard引入,随着Gaskell和Lau的简化对流界限准则(CBC)的引入而闻名,归一化变量图(NVD)是发展和分析HO和HR格式的有用工具。算法

NVF是基于局部归一化变量的面公式,用于构造面 f f 处的变量 ϕ f \phi_f ,该方法用到的是迎风( ϕ C \phi_C )、背风( ϕ D \phi_D )、和远迎风( ϕ U \phi_U )节点的值,以下图所示。架构

在这里插入图片描述

归一化变量的表达式为
ϕ ~ = ϕ ϕ U ϕ D ϕ U \tilde \phi=\frac{\phi-\phi_U}{\phi_D-\phi_U}
根据这个归一化关系, ϕ f \phi_f 的函数
ϕ f = f ( ϕ U , ϕ C , ϕ D ) \phi_f=f(\phi_U,\phi_C,\phi_D)
转变为
ϕ ~ f = f ( ϕ ~ C ) \tilde\phi_f=f(\tilde \phi_C)
之因此没有了 ϕ U \phi_U ϕ D \phi_D ,是由于他俩的归一化值变成了
ϕ ~ U = 0 ϕ ~ D = 1 \tilde \phi_U=0 \quad\quad \tilde \phi_D=1
那么 ϕ C \phi_C 的归一化值 ϕ ~ C \tilde\phi_C 变成了 ϕ \phi 场光滑性的衡量指标, ϕ ~ C \tilde\phi_C 的值在0到1之间( 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 )表明着单调廓线, ϕ ~ C \tilde\phi_C 的值小于0( ϕ ~ C < 0 \tilde\phi_C<0 )或者大于1( ϕ ~ C > 1 \tilde\phi_C>1 )则代表 C C 处存在极值。此外, ϕ ~ C 0 \tilde\phi_C\approx0 ϕ ~ C 1 \tilde\phi_C\approx1 代表是个梯度突跳(gradient jump)。这些构型都展现在下图中,图(a) ϕ ~ C < 0 \tilde\phi_C<0 ,极值;(b) ϕ ~ C > 1 \tilde\phi_C>1 ,极值;(c) ϕ ~ C = 1 \tilde\phi_C=1 ,突跳;(d) ϕ ~ C = 0 \tilde\phi_C=0 ,突跳;(e) 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 ,单调。
在这里插入图片描述app

在这里插入图片描述

归一化也可用于把高阶HO格式的函数关系转化为 ϕ ~ f \tilde\phi_f ϕ ~ C \tilde\phi_C 之间的线性关系,例如,前一章所讲到的HO格式的函数关系可转化成下面的归一化形式框架

格式 原形式 归一化形式
迎风Upwind ϕ f = ϕ C \phi_f=\phi_C ϕ ~ f = ϕ ~ C \Rightarrow\tilde\phi_f=\tilde\phi_C
中心差分CD ϕ f = 1 2 ( ϕ C + ϕ D ) \phi_f=\displaystyle\frac{1}{2}(\phi_C+\phi_D) ϕ ~ f = 1 2 ( 1 + ϕ ~ C ) \Rightarrow\tilde\phi_f=\displaystyle\frac{1}{2}(1+\tilde\phi_C)
二阶迎风SOU ϕ f = 3 2 ϕ C 1 2 ϕ U \phi_f=\displaystyle\frac{3}{2}\phi_C-\frac{1}{2}\phi_U ϕ ~ f = 3 2 ϕ ~ C \Rightarrow\tilde\phi_f=\displaystyle\frac{3}{2}\tilde\phi_C
FROMM ϕ f = ϕ C + 1 4 ( ϕ D ϕ U ) \phi_f=\phi_C+\displaystyle\frac{1}{4}(\phi_D-\phi_U) ϕ ~ f = ϕ ~ C + 1 4 \Rightarrow\tilde\phi_f=\tilde\phi_C+\displaystyle\frac{1}{4}
QUICK ϕ f = 3 4 ϕ C 1 8 ϕ U + 3 8 ϕ D \phi_f=\displaystyle\frac{3}{4}\phi_C-\frac{1}{8}\phi_U+\frac{3}{8}\phi_D ϕ ~ f = 3 8 + 3 4 ϕ ~ C \Rightarrow\tilde\phi_f=\displaystyle\frac{3}{8}+\displaystyle\frac{3}{4}\tilde\phi_C
背风Downwind ϕ f = ϕ D \phi_f=\phi_D ϕ ~ f = 1 \Rightarrow\tilde\phi_f=1

如此一来,对于基于三个节点值的全部高阶格式, ϕ ~ f \tilde\phi_f 总能表示成 ϕ ~ C \tilde\phi_C 的线性函数关系,即, ϕ ~ f = l ϕ ~ C + k \tilde\phi_f=l\tilde\phi_C+k ,其中 l l k k 是由格式所决定的系数。所以,若是把 ϕ ~ f \tilde\phi_f 做为 ϕ ~ C \tilde\phi_C 的函数绘制在 ( ϕ ~ C , ϕ ~ f ) (\tilde\phi_C,\tilde\phi_f) 平面上,那么这些格式的函数关系在图中将呈现为一条直线,所获得图形就是所谓的归一化变量图(Normalized Variable Diagram (NVD))。上述格式绘制出的NVD以下。ide

在这里插入图片描述

NVD代表,除了一阶迎风格式和背风格式,全部的二阶和三阶格式都经过了点 Q ( 0.5 , 0.75 ) Q(0.5,0.75) (均匀网格状况下)。实际上,若是一种格式具备二阶精度,那么它必需要经过 Q Q 点。此外,若是其在 Q Q 处的斜率是0.75,那么这种格式是三阶精度的(好比,QUICK格式)。迎风格式是很是扩散的,而背风格式是很是压缩的(反扩散)。所以,从NVD中能够形象地推出,若格式的函数关系与迎风格式的很是接近,那么它是扩散的,若是与背风格式的很是接近,那么它是压缩的。svg

前一章所讲的高阶格式比一阶迎风格式的截断偏差要小不少,还能保持稳定性。然而,这些格式的一个主要缺点是它们的无界性,即,它们倾向于在突跳和尖锐梯度处产生输运变量的过冲/下冲甚至是振荡。纵然在某些应用中,小的过冲或是振荡是尚可容忍的,然而,在某些情形下它们会致使灾难性的后果,好比在湍流计算中(输运变量多是湍流粘性,算不许就完蛋了)。函数

这些在尖锐梯度处的振荡特性是全部高阶线性对流格式的特性。实际上这些格式并非单调的,由于它们会产生局部最大或最小值,即,它们不能保持极值(extrema preserving)。若是一个格式能保持极值,那么解中的最大值必须不能再增大,且最小值不能再减少(该格式不能产生过冲/下冲)。实际上,Godunov和Ryabenki代表,任何单调的线性数值格式最多只能是一阶精度的,这也就意味着全部高阶线性格式不能保持单调性,而若要构造保持单调的格式,就须要使用非线性限制函数(non-linear limiter functions)。有了这样的理解,发展高阶无振荡对流格式的法子就能被分红两大组,第一类是,给一阶迎风格式添加上反扩散通量来加以限制(显然,是限制一阶迎风格式的假扩散特性),该添加项要使得产生的格式可以求解出尖锐的梯度而不会产生振荡才好;第二类是,给无界的高阶格式引入光滑扩散通量,来衰减非物理振荡。工具

因为它们的多步骤天性和在平衡两种通量时的困难,通量混合技术在数值算法中的代价是高昂的。这也就是为何在本书中,两种发展高分辨率算法的方法,只讲到通量限制方法。第一种是复合流程,其把高阶格式与有界低阶格式联合起来,且由一特定判据所控制的开关来切换两种格式。第二种方法是基于把扩散的一阶迎风项上添加一个反扩散通量,由通量限制器实现。此时,获得的高分辨率格式也称为总变差衰减(Total Variational Diminishing(TVD))格式,随后小节中将会详细讲述。

先讲解复合格式方法,用归一化变量公式(NVF)的架构来说解,并用归一化变量图(NVD)使讲解更具形象。所以,将首先讲解NVF和NVD,NVD的使用将做为一重要工具,来引出为确保高阶插值算法知足有界性,需定义什么准则。

2 对流有界准则(The Convection Boundness Criterion(CBC))

数值格式的预期效果应能保持其所描述或近似现象的物理特性,所以有界对流格式应该知足的条件,可经过分析对流的物理特性来更好理解。由于对流是把物理量从上游向下游输运的,那么对于对流的近似应拥有该输运特性。所以,数值对流格式应该是偏迎风的,否则的话它们缺少对流稳定性。那么,除了横跨界面两侧的节点 ϕ C \phi_C ϕ D \phi_D 以外,远迎风节点的值 ϕ U \phi_U 对于分析对流格式也是很是重要的,而再远的节点就意义较小了。在前面将的NVF中,值都作了归一化处理,也考虑了 ϕ U \phi_U 的影响,其很是好使,可有助于肯定数值对流格式在什么条件下是单调的。Spekreise和Barth和Jaspersen定义了一种单调格式(或者说是有界格式),该格式包含了面周围的全部邻居,Leonard和Gaskel和Lau基于他们的单调性定义,只给出了局部坐标系统上的邻近节点的定义形式,即
min ( ϕ C , ϕ D ) ϕ f max ( ϕ C , ϕ D ) \min(\phi_C,\phi_D)\le\phi_f\le\max(\phi_C,\phi_D)
归一化后,上式条件变为
min ( ϕ ~ C , 1 ) ϕ f max ( ϕ ~ C , 1 ) \min(\tilde\phi_C,1)\le\phi_f\le\max(\tilde\phi_C,1)
Gaskell和Lau提出的对于隐式稳态流动计算中的对流有界准则(The Convection Boundness Criterion(CBC))代表,一个拥有有界特性的格式,其函数关系应该是连续的,其界限应该是在 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 单调区间内位于下限 ϕ ~ C \tilde\phi_C 和上限 1 1 之间,应该经过点 ( 0 , 0 ) (0,0) ( 1 , 1 ) (1,1) ,应该在 ϕ ~ C > 1 \tilde\phi_C>1 ϕ ~ C < 0 \tilde\phi_C<0 区间的函数关系 f ( ϕ ~ C ) f(\tilde\phi_C) 等于 ϕ ~ C \tilde\phi_C ,上述条件在数学上的表达式为
ϕ ~ f = { f ( ϕ ~ C ) c o n t i n u o u s f ( ϕ ~ C ) = 1 i f   ϕ ~ C = 1 f ( ϕ ~ C )   w i t h   ϕ ~ C < f ( ϕ ~ C ) < 1 i f   0 < ϕ ~ C < 1 f ( ϕ ~ C ) = 0 i f   ϕ ~ C = 0 f ( ϕ ~ C ) = ϕ ~ C i f   ϕ ~ C < 0   o r   ϕ ~ C > 1 \tilde\phi_f=\begin{cases} f(\tilde\phi_C) & continuous \\ f(\tilde\phi_C)=1 & if~\tilde\phi_C=1 \\ f(\tilde\phi_C)~with~\tilde\phi_C<f(\tilde\phi_C)<1 & if~0<\tilde\phi_C<1 \\ f(\tilde\phi_C)=0 & if~\tilde\phi_C=0 \\ f(\tilde\phi_C)=\tilde\phi_C & if~\tilde\phi_C<0~or~\tilde\phi_C>1 \end{cases}
换句话说,知足上述关系的格式都是有界的格式,上述关系在NVD中的区域即下图的阴影部分。
在这里插入图片描述
在这里插入图片描述

结合上面俩图,就很是容易理解对流有界准则了。当 ϕ C \phi_C 是单调廓线的时候,在单元面处的插值廓线不该产生新的极值,因此其应该由面两侧的节点 ϕ \phi 值所限制,即 min ( ϕ C , ϕ D ) ϕ f max ( ϕ C , ϕ D ) \min(\phi_C,\phi_D)\le\phi_f\le\max(\phi_C,\phi_D) 。当 ϕ C \phi_C 逼近 ϕ D \phi_D 且仍在单调区间的时候, ϕ f \phi_f 也应该趋近于 ϕ D \phi_D 。当 ϕ C \phi_C 等于 ϕ D \phi_D 的时候, ϕ f \phi_f 也应该等于 ϕ D \phi_D ,此时 ( ϕ ~ C , ϕ ~ f ) (\tilde\phi_C,\tilde\phi_f) 经过点 ( 1 , 1 ) (1,1) 。当 ϕ C \phi_C 的值使得 ϕ ~ C > 1 \tilde\phi_C>1 ϕ f \phi_f 应给成迎风值,即, ϕ C \phi_C 。这样作的效果是,在知足了 ϕ f \phi_f 是由单元面两侧的节点值所限定的前提下,产生了最大可能的出流条件。如此一来,就会把错误的振荡衰减掉,由于 ϕ C \phi_C 将趋向于一个更低的值,得益于在这些条件中出流比入流要大的缘故(这句话的意思是说,对于单元 C C 来讲,此时其流出的量比流入的量要多,因此 ϕ C \phi_C 会变小,而不会出现振荡)。

所以,若是没有内在的物理机制(好比源项)来产生极值,那么极值将会消亡。一样的机理也对 ϕ ~ C < 0 \tilde\phi_C<0 的状况成立。随着在单调区域的 ϕ C \phi_C 逼近 ϕ U \phi_U ϕ f \phi_f 也会等于迎风值 ϕ C \phi_C 直到 ϕ C = ϕ U \phi_C=\phi_U ,代表 ϕ f \phi_f 的廓线要通过点 ( 0 , 0 ) (0,0) 。当 ϕ ~ C < 0 \tilde\phi_C<0 ϕ ~ C > 1 \tilde\phi_C>1 时,对流效应占优,因此这些区域选择迎风近似就恰如其分了。

3 高分辨率格式

使用NVD来构造有界高阶格式,即,高分辨率格式,是相对比较简单易行的。任何高阶格式可使用ad-hoc集合曲线(特殊曲线集,即分段直线所链接成的曲线,或者说是多点构造的曲线,两点间用直线表示)来令其具有有界性。

为了构造高分辨率格式,在区域 0 < ϕ ~ C < 1 0<\tilde\phi_C<1 的单调廓线应该经过点 ( 0 , 0 ) (0,0) ( 1 , 1 ) (1,1) ,同时须要在NVD中的上三角阴影区域。另外一方面,在非单调区域,即, ϕ ~ C < 0 \tilde\phi_C<0 ϕ ~ C > 1 \tilde\phi_C>1 区域,廓线应该使用迎风廓线。用这种方式构造的一些很是众所周知的高分辨率格式展现在下图中(关于MUSCL格式的图形和后面公式并不对应?鉴于书中所讲,应以公式为准)。

NVD of several HR schemes

为了提升收敛特性,任何复合高分辨率格式应该避免在其廓线链接点处出现较硬的角度,也不该该出现水平或是垂直的廓线。例如,对于SMART格式,它是由QUICK格式构造的,为了进一步提升其收敛性呢,能够作些微小的修正。即,把垂直廓线修正成区间 0 < ϕ ~ C < 1 / 6 0<\tilde\phi_C<1/6 ϕ ~ f = 3 ϕ ~ C \tilde\phi_f=3\tilde\phi_C ,对于STOIC格式,则是用 0 < ϕ ~ C < 1 / 5 0<\tilde\phi_C<1/5 区间。同时把SMART和STOIC格式中的水平廓线也做以修正,一种修正方法是把区间 7 / 10 < ϕ ~ C < 1 7/10<\tilde\phi_C<1 所对应的 9 / 10 < ϕ ~ f < 1 9/10<\tilde\phi_f<1 替换成 ϕ ~ f = ϕ ~ C / 3 + 2 / 3 \tilde\phi_f=\tilde\phi_C/3+2/3 。亦可用其它修正方法(好比,能够对复合廓线在 0.95 < ϕ ~ C < 1 0.95<\tilde\phi_C<1 的水平区间进行修正)。一样,也能够对有界CD格式作修正来提升其收敛特性。SMART、STOIC、SUPERBEE格式的修正NVD以下图所示。

在这里插入图片描述

数学上,这些复合高分辨率格式的函数关系为
MINMOD ϕ ~ f = { 3 2 ϕ ~ C 0 ϕ ~ C 1 2 1 2 ϕ ~ C + 1 2 1 2 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Bounded CD ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 ϕ ~ C e l s e w h e r e OSHER ϕ ~ f = { 3 2 ϕ ~ C 0 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e SMART ϕ ~ f = { 3 4 ϕ ~ C + 3 8 0 ϕ ~ C 5 6 1 5 6 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified SMART ϕ ~ f = { 3 ϕ ~ C 0 ϕ ~ C 1 6 3 4 ϕ ~ C + 3 8 1 6 ϕ ~ C 7 10 1 3 ϕ ~ C + 2 3 7 10 ϕ ~ C 1 ϕ ~ C e l s e w h e r e STOIC ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 2 3 4 ϕ ~ C + 3 8 1 2 ϕ ~ C 5 6 1 5 6 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified STOIC ϕ ~ f = { 3 ϕ ~ C 0 ϕ ~ C 1 5 1 2 ϕ ~ C + 1 2 1 5 ϕ ~ C 1 2 3 4 ϕ ~ C + 3 8 1 2 ϕ ~ C 7 10 1 3 ϕ ~ C + 2 3 7 10 ϕ ~ C 1 ϕ ~ C e l s e w h e r e MUSCL ϕ ~ f = { 2 ϕ ~ C 0 ϕ ~ C 1 4 ϕ ~ C + 1 4 1 4 ϕ ~ C 3 4 1 3 4 ϕ ~ C 1 ϕ ~ C e l s e w h e r e SUPERBEE ϕ ~ f = { 1 2 ϕ ~ C + 1 2 0 ϕ ~ C 1 2 3 2 ϕ ~ C 1 2 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e Modified SUPERBEE ϕ ~ f = { 2 ϕ ~ C 0 ϕ ~ C 1 3 1 2 ϕ ~ C + 1 2 1 3 ϕ ~ C 1 2 3 2 ϕ ~ C 1 2 ϕ ~ C 2 3 1 2 3 ϕ ~ C 1 ϕ ~ C e l s e w h e r e \begin{aligned} \text{MINMOD} \quad &\tilde\phi_f=\begin{cases} \frac{3}{2}\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{1}{2}\tilde\phi_C+\frac{1}{2} & \frac{1}{2}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Bounded CD} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2}\tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{OSHER} \quad &\tilde\phi_f=\begin{cases} \frac{3}{2}\tilde\phi_C & 0\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{SMART} \quad &\tilde\phi_f=\begin{cases} \frac{3}{4}\tilde\phi_C+\frac{3}{8} & 0\le\tilde\phi_C\le\frac{5}{6} \\ 1 & \frac{5}{6}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified SMART} \quad &\tilde\phi_f=\begin{cases} 3\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{6} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{6}\le\tilde\phi_C\le\frac{7}{10} \\ \frac{1}{3}\tilde\phi_C+\frac{2}{3} & \frac{7}{10}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{STOIC} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2}\tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{2}\le\tilde\phi_C\le\frac{5}{6} \\ 1 & \frac{5}{6}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified STOIC} \quad &\tilde\phi_f=\begin{cases} 3\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{5} \\ \frac{1}{2}\tilde\phi_C+\frac{1}{2} & \frac{1}{5}\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{4}\tilde\phi_C+\frac{3}{8} & \frac{1}{2}\le\tilde\phi_C\le\frac{7}{10} \\ \frac{1}{3}\tilde\phi_C+\frac{2}{3} & \frac{7}{10}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{MUSCL} \quad &\tilde\phi_f=\begin{cases} 2\tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{4} \\ \tilde\phi_C+\frac{1}{4} & \frac{1}{4}\le\tilde\phi_C\le\frac{3}{4} \\ 1 & \frac{3}{4}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{SUPERBEE} \quad &\tilde\phi_f=\begin{cases} \frac{1}{2} \tilde\phi_C+\frac{1}{2} & 0\le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{2} \tilde\phi_C & \frac{1}{2}\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \\ \text{Modified SUPERBEE} \quad &\tilde\phi_f=\begin{cases} 2 \tilde\phi_C & 0\le\tilde\phi_C\le\frac{1}{3} \\ \frac{1}{2} \tilde\phi_C+\frac{1}{2} & \frac{1}{3} \le\tilde\phi_C\le\frac{1}{2} \\ \frac{3}{2} \tilde\phi_C & \frac{1}{2}\le\tilde\phi_C\le\frac{2}{3} \\ 1 & \frac{2}{3}\le\tilde\phi_C\le1 \\ \tilde\phi_C & elsewhere \end{cases} \\ \end{aligned}
根据这种方法还发展了许多其它的格式,好比CLAM、UTOPIA、SHARP、ULTRA-SHARP等,不逐一列举。

4 TVD框架(总变差衰减框架)

另外一种发展高分辨率对流格式的流行方法是总变差衰减(Total Variation Diminishing(TVD))框架。在数值求解变量 ϕ \phi 的对流偏微分方程时,总变差(TV)定义为
T V = i ϕ i + 1 ϕ i TV=\sum_i|\phi_{i+1}-\phi_i|
其中 i i 表明空间求解域内的节点标识。若是解中的TV不随时间变化,那么该数值方法就被称为是总变差衰减的,数学上,等效于
T V ( ϕ t + Δ t ) T V ( ϕ t ) TV(\phi^{t+\Delta t}) \le TV(\phi^t)
在其影响深远的文章里,Harten证实了单调格式是TVD的,并且一个TVD格式是保持单调性的。一个保持单调性的格式并不会在求解域内产生任何新的局部极值,即,(原有的)局部最小值不会再减少,(原有的)局部最大值不会再增大。

这儿并不打算给出TVD方法的完整数学推导,而是,简单地讲解下构造TVD格式的方法,这里的方法是基于Sweby的工做。

回想下前面讲过的对流格式,一阶迎风格式具备强扩散性,而二阶中心差分格式具备高色散性。要是能有一种格式介于迎风和中心差分之间,即,拥有迎风格式的稳定性和中心差分格式的精度,岂不是很好。为构造这样的格式,从中心差分格式中寻找灵感
ϕ f = 1 2 ( ϕ D + ϕ C ) C D = ϕ C u p w i n d + 1 2 ( ϕ D ϕ C ) a n t i d i f f u s i v e   f l u x \phi_f=\underbrace{\frac{1}{2}(\phi_D+\phi_C)}_{CD}=\underbrace{\phi_C}_{upwind}+\underbrace{\frac{1}{2}(\phi_D-\phi_C)}_{anti-diffusive~flux}
这里的 C C D D 依旧表明着迎风节点和背风节点,而 f f 则表示节点 C C D D 之间的单元面。不难发现,中心差分格式能够写成迎风格式和一个通量的组合形式,该通量应该是起到了反扩散效应,从而把本来迎风格式的扩散性整成了中心差分格式的色散性。该通量的效果是使格式有了二阶精度,负面效应则是因为数值扩散效应的减弱而出现了非物理振荡。所以,一个较好的处理方式应该是,把这个反扩散通量的一部分添加到迎风格式中去,从而既能保持二阶精度,又不会产生非物理振荡。一种方法是把该通量乘上个限制函数(也叫限制器或通量限制器),在振荡可能发生的区域(如,较大梯度的地方)防止其过分使用(好比直接把它设成0),而在光滑区域发挥其最大贡献。把这种限制器定义为 ψ ( r ) \psi(r) ,其中 r r 常取为两个连续梯度的比值,这样, ϕ f \phi_f 的计算为
ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) w i t h   r f = ϕ C ϕ U ϕ D ϕ C \phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)\quad\quad with~r_f=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
其中 U U 是远迎风点,而 D D 是背风点, C C 仍为迎风点。为了保证反扩散通量的符号特性, ψ ( r f ) \psi(r_f) 须要是非负的,即 ψ ( r f ) 0 \psi(r_f)\ge0 (否则起不到反扩散效果啊)。此外,若是 r f < 0 r_f<0 则代表节点 C C 是局部极值,而 r f = 0 r_f=0 代表有局部的陡峭梯度(有个台阶,突跳状况),此时应该让 ψ ( r f ) = 0 \psi(r_f)=0 ,这样才能不让反扩散通量起做用,防止出现非物理振荡。那么,发展TVD格式就变成了找到限制器 ψ ( r f ) \psi(r_f) ,使得数值格式是TVD或者单调的。这些限制器必需要知足哪些条件,才能保证对流格式是保单调的呢?

感受书中这部分(P444-446)的理论推导过程并不正确,但结论毋庸置疑是没错的,因为正确的推导我也没搞清楚,因此直接把结论给出来吧。那就是这个限制器 ψ ( r ) \psi(r) 必需要知足 r > 0 r>0 时(为方便计,把下标 f f 去掉了)
ψ ( r ) 2 a n d ψ ( r ) 2 r \psi(r)\le2\quad and \quad \psi(r)\le2r
与前面提到的 r 0 r\le0 ψ ( r ) = 0 \psi(r)=0 联合起来,给出TVD格式应该知足的条件,即限制器应该知足的条件为
ψ ( r ) = { m i n ( 2 r , 2 ) r > 0 0 r 0 \psi(r)=\begin{cases} min(2r,2) & r>0 \\ 0 & r\le0 \end{cases}
在这里插入图片描述

如上图所示,这个条件能够画在 r ψ r-\psi 图上,这也被称为Sweby图,图中的蓝色区域即是TVD单调区域。使用这个图的话,就很容易掌握TVD格式的公式。任何通量限制器的公式只要是能落在TVD单调区域内的都能产生出TVD格式来,Sweby图跟前面讲的NVD是很是相似的。

前面讲过的全部格式的限制器也能够推导出来,还能够把它们的函数关系绘制在Sweby图中。显然,中心差分格式的限制器为 ψ C D = 1 \psi^{CD}=1 ,而二阶迎风格式的限制器为
ϕ f = ϕ C + 1 2 ψ S O U ( r f ) ( ϕ D ϕ C ) = 3 2 ϕ C 1 2 ϕ U ψ S O U ( r f ) = ϕ C ϕ U ϕ D ϕ C = r f \begin{aligned} \phi_f&=\phi_C+\frac{1}{2}\psi^{SOU}(r_f)(\phi_D-\phi_C)=\frac32\phi_C-\frac12\phi_U \\ \Rightarrow \psi^{SOU}(r_f)&=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}=r_f \end{aligned}
中心差分格式和二阶迎风格式的限制器展现在Sweby图中,以下,CD为中心差分格式,SOU为二阶迎风格式。
在这里插入图片描述

Sweby还注意到,因为 r 0 r\le0 ψ ( r ) = 0 \psi(r)=0 ,因此在解中有极值的地方就没有二阶精度了。二阶迎风和中心差分格式都是二阶格式,从图中可看到,他俩的限制器都通过了点 ( 1 , 1 ) (1,1) 。此外,Van Leer的工做代表,任何二阶格式都能写成中心差分格式和二阶迎风格式的加权平均形式。因此呢,若是一个格式有二阶精度,那么它的限制器必需要穿过点 ( 1 , 1 ) (1,1) 才对,并且,该限制器应该位于中心差分格式和二阶迎风格式限制器所围起来的区域,即上图中的蓝色区域才对。若是要用上一小节的NVD(归一化变量图)来表示该区域,那么是下面这样子的。
在这里插入图片描述

一样也能够求出来其它高阶格式的限制器来,即
{ U p w i n d ψ ( r f ) = 0 D o w n w i n d ψ ( r f ) = 2 F R O M M ψ ( r f ) = 1 + r f 2 S O U ψ ( r f ) = r f C D ψ ( r f ) = 1 Q U I C K ψ ( r f ) = 3 + r f 4 \begin{cases} Upwind & \psi(r_f)=0 \\ Downwind & \psi(r_f)=2 \\ FROMM & \psi(r_f)=\dfrac{1+r_f}{2} \\ SOU & \psi(r_f)=r_f \\ CD & \psi(r_f)=1 \\ QUICK & \psi(r_f)=\dfrac{3+r_f}{4} \end{cases}
这些限制器的函数关系以下图所示,除了迎风格式的限制器外,全部其它格式的限制器都并不是彻底位于单调区间内,这也就是为何这些格式都是无界的。
在这里插入图片描述

经过把这些不一样格式的限制器 ψ ( r f ) \psi(r_f) 限制到单调区间内,即可把这些高阶格式转化为高分辨率TVD格式。经过这种方式能够发展不少种TVD格式,它们的限制器以下图所示,它们限制器的函数关系以下。

在这里插入图片描述

{ S U P E R B E E ψ ( r f ) = max ( 0 , min ( 1 , 2 r f ) , min ( 2 , r f ) ) M I N M O D ψ ( r f ) = max ( 0 , min ( 1 , r f ) ) O S H E R ψ ( r f ) = max ( 0 , min ( 2 , r f ) ) V a n   L e e r ψ ( r f ) = r f + r f 1 + r f M U S C L ψ ( r f ) = max ( 0 , min ( 2 r f , ( r f + 1 ) / 2 , 2 ) ) \begin{cases} SUPERBEE & \psi(r_f)= \max(0,\min(1,2r_f),\min(2,r_f)) \\ MINMOD & \psi(r_f)=\max(0,\min(1,r_f)) \\ OSHER & \psi(r_f)= \max(0,\min(2,r_f)) \\ Van~Leer & \psi(r_f)=\dfrac{r_f+|r_f|}{1+|r_f|} \\ MUSCL & \psi(r_f)= \max(0,\min(2r_f,(r_f+1)/2, 2)) \end{cases}

5 NVF-TVD关系

既然NVF和TVD公式都是用来添加有界性的,只不过是采用不一样的方法而已,那么它俩应该是有某种关联才对啊。下面先推导出 r f r_f ϕ ~ C \tilde\phi_C 的关系式,而后比较NVF-CBC和TVD-CBC,最后导出通用的转化关系,可让任何TVD格式的函数关系写成NVF框架的形式,或者是反过来,把NVF格式写成TVD形式。

r f r_f ϕ ~ C \tilde\phi_C 的关系不难推导,直接由它们的定义可得
r f = ϕ C ϕ U ϕ D ϕ C ϕ ~ C = ϕ C ϕ U ϕ D ϕ U } r f = ( ϕ C ϕ U ) / ( ϕ D ϕ U ) ( ϕ D ϕ U + ϕ U ϕ C ) / ( ϕ D ϕ U ) = ϕ ~ C 1 ϕ ~ C ϕ ~ C = r f 1 + r f \left. \begin{aligned} r_f&=\frac{\phi_C-\phi_U}{\phi_D-\phi_C} \\ \tilde\phi_C&=\frac{\phi_C-\phi_U}{\phi_D-\phi_U} \end{aligned} \right\}\Rightarrow r_f=\frac{(\phi_C-\phi_U)/(\phi_D-\phi_U)}{(\phi_D-\phi_U+\phi_U-\phi_C)/(\phi_D-\phi_U)}= \frac{\tilde\phi_C}{1-\tilde\phi_C}\\ \Rightarrow \tilde\phi_C=\frac{r_f}{1+r_f}
使用这个式子,能够比较两种框架下的线性格式,好比迎风格式在TVD框架下的限制器为 ψ ( r f ) = 0 \psi(r_f)=0 ,那么其应该等效于NVF框架下的迎风格式 ϕ ~ f = ϕ ~ C \tilde\phi_f=\tilde\phi_C ,二者的等效性很容易证实
ψ ( r f ) = 0 ϕ f = ϕ C ϕ ~ f = ϕ ~ C \psi(r_f)=0\Rightarrow\phi_f=\phi_C\Rightarrow\tilde\phi_f=\tilde\phi_C
此外,无论是何种格式,在 r f 0 r_f\le0 时,TVD-CBC中的限制器 ψ ( r f ) = 0 \psi(r_f)=0 变为纯粹的迎风格式,有
r f 0 ϕ ~ C 1 ϕ ~ C 0 { ϕ ~ C 0 ϕ ~ C 1 r_f\le0 \Rightarrow \frac{\tilde\phi_C}{1-\tilde\phi_C}\le0 \Rightarrow \left\{\begin{aligned} \tilde\phi_C\le0 \\ \tilde\phi_C\ge1 \end{aligned}\right.
这跟在NVF-CBC中的迎风格式的使用区间是彻底一致的,即 ϕ ~ C 0 \tilde\phi_C\le0 ϕ ~ C 1 \tilde\phi_C\ge1 时,使用纯粹的迎风格式。

还有,在NVF-CBC中,函数关系在区间 0 ϕ ~ C 1 0\le\tilde\phi_C\le1 是单调增长的,而在Sweby图中该区域则是 0 r f + 0\le r_f \le +\infin ,这两个区间其实是一样的,以下所示:
ϕ ~ C 1 r f = ϕ ~ C 1 ϕ ~ C + \tilde\phi_C \rightarrow 1 \Rightarrow r_f=\frac{\tilde\phi_C}{1-\tilde\phi_C} \rightarrow +\infin
更进一步,对于TVD-CBC的条件
ψ ( r f ) 2 \psi(r_f)\le2
等效的NVF-CBC条件为
ψ ( r f ) = 2 ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) } ϕ f = ϕ C + ( ϕ D ϕ C ) = ϕ D ϕ ~ f = 1 \left. \begin{aligned} &\psi(r_f)=2\\ &\phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C) \end{aligned} \right\} \Rightarrow \phi_f=\phi_C+(\phi_D-\phi_C)=\phi_D \Rightarrow \tilde\phi_f=1

ψ ( r f ) 2 ϕ ~ f 1 \psi(r_f)\le2 \Rightarrow \tilde\phi_f\le1
正是NVF-CBC所应知足的条件。TVD-CBC的限制器 ψ ( r f ) \psi(r_f) 应知足的最后一个条件是
ψ ( r f ) 2 r f \psi(r_f) \le 2r_f
等效的NVF-CBC条件为
ψ ( r f ) = 2 r f ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) } ϕ f = ϕ C + ϕ C ϕ U ϕ D ϕ C ( ϕ D ϕ C ) = 2 ϕ C ϕ U \left. \begin{aligned} & \psi(r_f) = 2r_f\\ &\phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C) \end{aligned} \right\} \Rightarrow \phi_f=\phi_C+\frac{\phi_C-\phi_U}{\phi_D-\phi_C}(\phi_D-\phi_C)=2\phi_C-\phi_U
归一化后,为
ϕ f = 2 ϕ C ϕ U ϕ f ϕ U = 2 ϕ C 2 ϕ U ϕ ~ f = 2 ϕ ~ C \phi_f=2\phi_C-\phi_U \Rightarrow \phi_f-\phi_U=2\phi_C-2\phi_U \Rightarrow \tilde\phi_f=2\tilde\phi_C
这比NVF-CBC要更加严格,这是两种方法的惟一不一样之处。基于该条件,TVD-CBC和修正的NVF-CBC(用TVD条件表示的NVF)以下图a所示,修正NVF-CBC的单调区间减少成了迎风线和蓝色线所围的区域;而修正TVD-CBC(用NVF条件表示的TVD)和NVF-CBC以下图b所示(NVF-CBC中的线 ϕ ~ C = 0 \tilde\phi_C=0 对应TVD-CBC中的线 r f = 0 r_f=0 )。再考虑到二阶精度的要求,TVD格式中的二阶精度要求必须经过点 ( 1 , 1 ) (1,1) ,即, ψ ( 1 ) = 1 \psi(1)=1 ,那么其等效的NVF值为
r f = 1 ϕ ~ C 1 ϕ ~ C = 1 ϕ ~ C = 0.5 ψ ( 1 ) = 1 ϕ f = ϕ C + 1 2 ψ ( 1 ) ( ϕ D ϕ C ) = 1 2 ( ϕ D + ϕ C ) ϕ f ϕ U = 1 2 ( ϕ D ϕ U + ϕ C ϕ U ) ϕ ~ f = 1 2 ( 1 + ϕ ~ C ) = 0.75 \begin{aligned} r_f=1 &\Rightarrow \frac{\tilde\phi_C}{1-\tilde\phi_C}=1 \Rightarrow \tilde\phi_C=0.5 \\ \psi(1)=1 &\Rightarrow \phi_f=\phi_C+\frac{1}{2}\psi(1)(\phi_D-\phi_C)= \frac{1}{2}(\phi_D+\phi_C) \\ & \Rightarrow \phi_f-\phi_U=\frac{1}{2}(\phi_D-\phi_U+\phi_C-\phi_U) \Rightarrow \tilde\phi_f=\frac{1}{2}(1+\tilde\phi_C)=0.75 \end{aligned}
这偏偏就是NVF中的点 ( 0.5 , 0.75 ) (0.5,0.75) 。前面提到过,Van Leer证实了任何二阶格式均可以写成中心差分格式和二阶迎风格式的加权平均形式,所以其函数关系应该位于CD和SOU格式的函数关系之间,以下图c所示。

至此,充分说明了NVF-CBC和TVD-CBC实际上就是一回事,即,异曲同工。
在这里插入图片描述

上面的过程也能够用来把TVD格式转化为等效的NVF格式,或者是反过来作转化。用NVF框架下的格式为例,面 f f 上的值 ϕ f \phi_f 表示为
ϕ f = f ( ϕ ~ C ) ( ϕ D ϕ U ) + ϕ U ϕ ~ C = ϕ C ϕ U ϕ D ϕ C \phi_f=f(\tilde\phi_C)(\phi_D-\phi_U)+\phi_U\quad\quad\tilde\phi_C=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
TVD格式的 ϕ f \phi_f 表达式为
ϕ f = ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) r f = ϕ C ϕ U ϕ D ϕ C \phi_f=\phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)\quad\quad r_f=\frac{\phi_C-\phi_U}{\phi_D-\phi_C}
结合上面两式,得
ϕ C + 1 2 ψ ( r f ) ( ϕ D ϕ C ) = f ( ϕ ~ C ) ( ϕ D ϕ U ) + ϕ U \phi_C+\frac{1}{2}\psi(r_f)(\phi_D-\phi_C)=f(\tilde\phi_C)(\phi_D-\phi_U)+\phi_U

ψ ( r f ) ( ϕ D ϕ C ) ( ϕ D ϕ U ) = 2 f ( ϕ ~ C ) ( ϕ D ϕ U ) ( ϕ D ϕ U ) 2 ( ϕ C ϕ U ) ( ϕ D ϕ U ) = 2 ( f ( ϕ ~ C ) ϕ ~ C ) \psi(r_f)\frac{(\phi_D-\phi_C)}{(\phi_D-\phi_U)}= 2f(\tilde\phi_C)\frac{(\phi_D-\phi_U)}{(\phi_D-\phi_U)}- 2\frac{(\phi_C-\phi_U)}{(\phi_D-\phi_U)} =2(f(\tilde\phi_C)-\tilde\phi_C)
上式左端项能够写为
ψ ( r f ) ( ϕ D ϕ C ) ( ϕ D ϕ U ) = ψ ( r f ) ( ϕ D ϕ U ϕ C + ϕ U ) ( ϕ D ϕ U ) = ψ ( r f ) ( 1 ϕ ~ C ) \psi(r_f)\frac{(\phi_D-\phi_C)}{(\phi_D-\phi_U)}= \psi(r_f)\frac{(\phi_D-\phi_U-\phi_C+\phi_U)}{(\phi_D-\phi_U)}= \psi(r_f)(1-\tilde\phi_C)
因而
ψ ( r f ) ( 1 ϕ ~ C ) = 2 ( f ( ϕ ~ C ) ϕ ~ C ) ψ ( r f ) = 2 f ( ϕ ~ C ) ϕ ~ C 1 ϕ ~ C \psi(r_f)(1-\tilde\phi_C)=2(f(\tilde\phi_C)-\tilde\phi_C) \Rightarrow \psi(r_f)=2\frac{f(\tilde\phi_C)-\tilde\phi_C}{1-\tilde\phi_C}
上式也可写成
f ( ϕ ~ C ) = ψ ( r f ) + 2 r f 2 ( 1 + r f ) f(\tilde\phi_C)=\frac{\psi(r_f)+2r_f}{2(1+r_f)}
这样便推得了NVF格式与TVD格式的转换关系。

举个例子,迎风格式在NVF框架下的函数关系是 ϕ ~ f = ϕ ~ C \tilde\phi_f=\tilde\phi_C ,则其TVD限制器为
ϕ ~ f = ϕ ~ C ψ ( r f ) = 2 f ( ϕ ~ C ) ϕ ~ C 1 ϕ ~ C = 2 ϕ ~ C ϕ ~ C 1 ϕ ~ C = 0 \tilde\phi_f=\tilde\phi_C \Rightarrow \psi(r_f)=2\frac{f(\tilde\phi_C)-\tilde\phi_C}{1-\tilde\phi_C}= 2\frac{\tilde\phi_C-\tilde\phi_C}{1-\tilde\phi_C}=0