有限体积法（10）——格式精度与待定系数法 CFD

前面用到的格式及其计算精度：

格式名称	公式	精度
中心差分格式	( ? ? ? x ) e = ? E ? ? P Δ x \left (\frac{\partial \phi}{\partial x}\right)_e=\frac{\phi_{E}-\phi_{P}}{\Delta x} (?x???)e?=Δx?E???P??	二阶
迎风格式	? w = ? W , ? e = ? P \phi_w=\phi_{W},\phi_e=\phi_P ?w?=?W?,?e?=?P?（流动沿 x x x正方向时）	一阶
QUICK格式	? e = 3 8 ? E + 6 8 ? P ? 1 8 ? W \phi_{e}=\frac{3}{8}\phi_{E}+\frac{6}{8}\phi_P-\frac{1}{8}\phi_{W} ?e?=83??E?+86??P??81??W?	三阶

下面将推导上述格式的误差和精度，以及介绍使用待定系数法构造差分格式的过程。
格式精度

文章图片

以一维情况为例，函数 ? ( x ) \phi(x) ?(x)在节点i i i 附近的泰勒级数展开为
? ( x + Δ x ) = ? ( x ) + ( ? ? ? x ) x Δ x + ( ? 2 ? ? x 2 ) x Δ x 2 2 + . . . (1) \begin{aligned} \phi(x+\Delta x)=\phi(x)+\left( \frac{\partial \phi}{\partial x} \right)_x \Delta x+\left( \frac{\partial^2\phi}{\partial x^2} \right)_x \frac{\Delta x^2}{2}+... \end{aligned} \tag{1} ?(x+Δx)=?(x)+(?x???)x?Δx+(?x2?2??)x?2Δx2?+...?(1)
使用上图中的标记以及各节点对应关系，则泰勒展开式写为
? E = ? P + ( ? ? ? x ) P Δ x + ( ? 2 ? ? x 2 ) P Δ x 2 2 + . . . (2) \begin{aligned} \phi_E=\phi_P+\left( \frac{\partial \phi}{\partial x} \right)_P \Delta x+\left( \frac{\partial^2\phi}{\partial x^2} \right)_P \frac{\Delta x^2}{2}+... \end{aligned} \tag{2} ?E?=?P?+(?x???)P?Δx+(?x2?2??)P?2Δx2?+...?(2)
整理上式有
( ? ? ? x ) P Δ x = ? E ? ? P ? ( ? 2 ? ? x 2 ) P Δ x 2 2 ? . . . ( ? ? ? x ) P = ? E ? ? P Δ x ? ( ? 2 ? ? x 2 ) P Δ x 2 ? . . . (3) \begin{aligned} &\left( \frac{\partial \phi}{\partial x} \right)_P \Delta x=\phi_E - \phi_P - \left( \frac{\partial^2 \phi}{\partial x^2}\right)_P \frac{\Delta x^2}{2}-... \\\\ &\left( \frac{\partial \phi}{\partial x} \right)_P=\frac{\phi_E-\phi_P}{\Delta x} -\left( \frac{\partial^2 \phi}{\partial x^2}\right)_P \frac{\Delta x}{2}-... \end{aligned} \tag{3} ?(?x???)P?Δx=?E???P??(?x2?2??)P?2Δx2??...(?x???)P?=Δx?E???P???(?x2?2??)P?2Δx??...?(3)
所以有
( ? ? ? x ) P = ? E ? ? P Δ x + 截断项 (4) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_P=\frac{\phi_E-\phi_P}{\Delta x} + 截断项 \end{aligned} \tag{4} (?x???)P?=Δx?E???P??+截断项?(4)
从公式 ( 3 ) (3) (3)可以看出，截断项是某个数与 Δ x \Delta x Δx的乘积，如果略去截断项，有
( ? ? ? x ) P ≈ ? E ? ? P Δ x (5) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_P \approx \frac{\phi_E-\phi_P}{\Delta x} \end{aligned} \tag{5} (?x???)P?≈Δx?E???P???(5)
使用公式 ( 5 ) (5) (5)来近似导数包含了一定的误差，该误差就是省略了截断项所导致的。根据公式 ( 3 ) (3) (3)可以看出，通过减小 Δ x \Delta x Δx的量来降低公式 ( 5 ) (5) (5)截断误差。一般而言，有限差分格式的截断项都含有 Δ x n \Delta x^n Δxn， Δ x \Delta x Δx的指数就叫做差分近似的阶数，它决定了当网格加密（ Δ x \Delta x Δx变小）时，截断误差趋于零的速率。因此，公式 ( 5 ) (5) (5)所表示的差分格式就是一阶精度的格式，我们写为
( ? ? ? x ) P = ? E ? ? P Δ x + O ( Δ x ) (6) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_P=\frac{\phi_E-\phi_P}{\Delta x} + O(\Delta x) \end{aligned} \tag{6} (?x???)P?=Δx?E???P??+O(Δx)?(6)
因为上式使用了节点 E E E和 P P P来计算节点 P P P处的导数 ( ? ? ? x ) P \left(\frac{\partial \phi}{\partial x}\right)_P (?x???)P?，公式 ( 6 ) (6) (6)就叫做前向差分格式。同理，我们可以推导出节点 P P P处的后向差分，
? ( x ? Δ x ) = ? ( x ) ? ( ? ? ? x ) x Δ x + ( ? 2 ? ? x 2 ) x Δ x 2 2 + . . . (7) \begin{aligned} \phi(x-\Delta x)=\phi(x)-\left( \frac{\partial \phi}{\partial x} \right)_x \Delta x+\left( \frac{\partial^2\phi}{\partial x^2} \right)_x \frac{\Delta x^2}{2}+... \end{aligned} \tag{7} ?(x?Δx)=?(x)?(?x???)x?Δx+(?x2?2??)x?2Δx2?+...?(7)
节点 P P P处的后向差分格式为
( ? ? ? x ) P = ? P ? ? W Δ x + O ( Δ x ) (8) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_P=\frac{\phi_P-\phi_W}{\Delta x} + O(\Delta x) \end{aligned} \tag{8} (?x???)P?=Δx?P???W??+O(Δx)?(8)
与前向差分相同，后向差分也是一阶精度的。前向和后向差分格式都只使用了两个节点值。
中心差分格式用公式 ( 1 ) (1) (1)减去公式 ( 7 ) (7) (7)得
? ( x + Δ x ) ? ? ( x ? Δ x ) = 2 ( ? ? ? x ) P Δ x + ( ? 3 ? ? x 3 ) P Δ x 3 3 ! + . . . (9) \begin{aligned} \phi(x+\Delta x)-\phi(x-\Delta x)= 2 \left( \frac{\partial \phi}{\partial x} \right)_P \Delta x + \left( \frac{\partial^3 \phi}{\partial x^3} \right)_P \frac{\Delta x^3}{3!} +... \end{aligned} \tag{9} ?(x+Δx)??(x?Δx)=2(?x???)P?Δx+(?x3?3??)P?3!Δx3?+...?(9)
整理之，然后得到 P P P点导数的另一个差分格式
( ? ? ? x ) P = ? E ? ? W 2 Δ x + O ( Δ x 2 ) (10) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_P = \frac{\phi_E-\phi_W}{2\Delta x} + O(\Delta x^2) \end{aligned} \tag{10} (?x???)P?=2Δx?E???W??+O(Δx2)?(10)
公式 ( 10 ) (10) (10)使用了节点 E E E和 W W W的值来计算中间节点 P P P处的导数，该差分格式称为中心差分格式。由其截断项可见，中心差分格式是二阶精度的。当网格尺度减小时，省略截断项导致的误差会以二次方的速率趋于零，这比一阶精度要快。
在前面的方程离散推导中，导数的计算主要存在与扩散项，并且求解导数的位置一般是在单元的界面处，例如 e e e或 w w w，那时用的中心差分格式为
( ? ? ? x ) e = ? E ? ? P Δ x = ? E ? ? P 2 ( Δ x / 2 ) (11) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_e = \frac{\phi_E-\phi_P}{\Delta x} =\frac{\phi_E-\phi_P}{2(\Delta x/2)} \end{aligned} \tag{11} (?x???)e?=Δx?E???P??=2(Δx/2)?E???P???(11)
根据上图可知，公式 ( 11 ) (11) (11)就是使用了节点 E E E和 P P P的值计算了其中点 e e e处的导数，对于均匀网格来说，它就是具有二阶精度的中心差分格式。
迎风格式类似公式 ( 4 ) 、 ( 5 ) 、 ( 6 ) (4)、(5)、(6) (4)、(5)、(6)，公式 ( 2 ) (2) (2)可以写成
? E = ? P + 截断项 (12) \begin{aligned} \phi_E = \phi_P + 截断项 \end{aligned} \tag{12} ?E?=?P?+截断项?(12)
如果省略截断项，用节点 P P P的值来近似表示节点 E E E的值，则就得出迎风格式
? E = ? P + O ( Δ x ) (13) \begin{aligned} \phi_E = \phi_P + O(\Delta x) \end{aligned} \tag{13} ?E?=?P?+O(Δx)?(13)
在对流项的离散中，需要近似的是单元界面上的变量值，则对于均匀网格有
? e = ? P + O ( Δ x / 2 ) = ? P + O ( Δ x ) ? w = ? W + O ( Δ x / 2 ) = ? P + O ( Δ x ) (14) \begin{aligned} \phi_e = \phi_P + O(\Delta x/2) = \phi_P + O(\Delta x) \\\\ \phi_w = \phi_W + O(\Delta x/2) = \phi_P + O(\Delta x) \end{aligned} \tag{14} ?e?=?P?+O(Δx/2)=?P?+O(Δx)?w?=?W?+O(Δx/2)=?P?+O(Δx)?(14)
【有限体积法（10）——格式精度与待定系数法】可见，迎风格式是一阶精度的。（这里说的迎风格式指的是基本迎风格式，即用上游值代替边界值，迎风格式也有其他高阶格式）
QUICK格式以用QUICK格式计算 ? e \phi_e ?e?为例，近似公式为
? e = 3 8 ? E + 6 8 ? P ? 1 8 ? W (15) \begin{aligned} \phi_{e}=\frac{3}{8}\phi_{E}+\frac{6}{8}\phi_P-\frac{1}{8}\phi_{W} \end{aligned} \tag{15} ?e?=83??E?+86??P??81??W??(15)
这里点 E 、 P 、 W E、P、W E、P、W称为模板点，我们希望通过模板点的值来近似界面e处的值，那么可以把近似公式写成待定系数的形式
? e = ( a? E + b? P ? c? W ) + O ( Δ x k ) (16) \begin{aligned} \phi_{e} = (a\ \phi_{E}+b\ \phi_P-c\ \phi_{W}) + O(\Delta x^k) \end{aligned} \tag{16} ?e?=(a ?E?+b ?P??c ?W?)+O(Δxk)?(16)
其中， a 、 b 、 c a、b、c a、b、c为待定系数。
函数 ? \phi ?在界面 e e e处的泰勒级数展开有
? E = ? e + Δ x 2 ( ? ? ? x ) e + 1 2 ( Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e + O ( Δ x 3 ) ? P = ? e ? Δ x 2 ( ? ? ? x ) e + 1 2 ( ? Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e + O ( Δ x 3 ) ? W = ? e ? 3 Δ x 2 ( ? ? ? x ) e + 1 2 ( ? 3 Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e + O ( Δ x 3 ) (17) \begin{aligned} \phi_E=\phi_e+ \frac{\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_e+ \frac{1}{2} \left ( \frac{\Delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_e +O(\Delta x^3) \\ \\ \phi_P=\phi_e- \frac{\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_e+ \frac{1}{2} \left ( -\frac{\Delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_e +O(\Delta x^3) \\ \\ \phi_W=\phi_e- \frac{3\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_e+ \frac{1}{2} \left ( -\frac{3\Delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_e +O(\Delta x^3) \\ \\ \end{aligned} \tag{17} ?E?=?e?+2Δx?(?x???)e?+21?(2Δx?)2(?x2?2??)e?+O(Δx3)?P?=?e??2Δx?(?x???)e?+21?(?2Δx?)2(?x2?2??)e?+O(Δx3)?W?=?e??23Δx?(?x???)e?+21?(?23Δx?)2(?x2?2??)e?+O(Δx3)?(17)
把公式 ( 17 ) (17) (17)代入到公式 ( 16 ) (16) (16)得
? e = a[ ? e + Δ x 2 ( ? ? ? x ) e + 1 2 ( Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e ] + b[ ? e ? Δ x 2 ( ? ? ? x ) e + 1 2 ( ? Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e ] + c[ ? e ? 3 Δ x 2 ( ? ? ? x ) e + 1 2 ( ? 3 Δ x 2 ) 2 ( ? 2 ? ? x 2 ) e ] + O ( Δ x 3 ) (18) \begin{aligned} \phi_e = &a\ \left [ \phi_e + \frac{\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_e + \frac{1}{2} \left(\frac{\Delta x}{2} \right)^2 \left( \frac{\partial^2 \phi}{\partial x^2} \right)_e \right] \\\\ &+b\ \left [ \phi_e - \frac{\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right) _e+ \frac{1}{2} \left(-\frac{\Delta x}{2} \right)^2 \left( \frac{\partial^2 \phi}{\partial x^2} \right)_e \right] \\\\ &+c\ \left [ \phi_e - \frac{3\Delta x}{2} \left( \frac{\partial \phi}{\partial x} \right) _e+ \frac{1}{2} \left(-\frac{3\Delta x}{2} \right)^2 \left( \frac{\partial^2 \phi}{\partial x^2} \right)_e \right] \\\\ &+ O(\Delta x^3) \end{aligned} \tag{18} ?e?=?a [?e?+2Δx?(?x???)e?+21?(2Δx?)2(?x2?2??)e?]+b [?e??2Δx?(?x???)e?+21?(?2Δx?)2(?x2?2??)e?]+c [?e??23Δx?(?x???)e?+21?(?23Δx?)2(?x2?2??)e?]+O(Δx3)?(18)
因为函数 ? \phi ?是不确定的，将函数 ? \phi ?及其各阶导数看作是变量，根据等号两边对应系数相等有
{ a + b + c = 1 1 2 a ? 1 2 b ? 3 2 c = 0 1 8 a + 1 8 b + 9 8 c = 0 (19) \left \{ \begin{aligned} a+b+c=1\\\\ \frac{1}{2}a-\frac{1}{2}b-\frac{3}{2}c=0\\\\ \frac{1}{8}a+\frac{1}{8}b+\frac{9}{8}c=0 \end{aligned} \right. \tag{19} ??????????????????????a+b+c=121?a?21?b?23?c=081?a+81?b+89?c=0?(19)
解方程组得
{ a = 3 8 b = 6 8 b = ? 1 8 (20) \left \{ \begin{aligned} a&=\frac{3}{8}\\\\ b&=\frac{6}{8}\\\\ b&=-\frac{1}{8} \end{aligned} \right. \tag{20} ????????????????????????abb?=83?=86?=?81??(20)
∴ \therefore ∴
? e = ( 3 8 ? E + 6 8 ? P ? 1 8 ? W ) + O ( Δ x 3 ) (21) \begin{aligned} \phi_{e} = ( \frac{3}{8}\phi_{E}+\frac{6}{8} \phi_P-\frac{1}{8} \phi_{W}) + O(\Delta x^3) \end{aligned} \tag{21} ?e?=(83??E?+86??P??81??W?)+O(Δx3)?(21)
可见QUICK格式具有三阶精度。上述过程就是使用待定系数法构建差分格式的一般过程，这里构建的是近似函数值的QUICK格式，也可以用于构建近似导数的差分格式。公式 ( 19 ) (19) (19)中方程式是根据 ? \phi ?的零阶、一阶和二阶导数的系数在等号两边对应相等得出的，其实如果把泰勒级数展开到三阶导数项也可以得出第四个方程，这样方程组就成了超定方程组了，这种情况方程组一般无解，也就是说，我们不能用 ? E 、 ? P 、 ? W \phi_E、\phi_P、\phi_W ?E?、?P?、?W?三个节点值构造出 ? e \phi_e ?e?的四阶差分近似。一般情况下，对于均匀网格，用 k k k个连续的模板点来近似 ? e \phi_e ?e?的精度最高为 k k k阶。（这里仅指对函数值的近似，近似导数时阶数会更低，例如下面用三个点近似一阶导数就是二阶精度的。）
用待定系数法构造导数的近似格式

文章图片

在QUICK格式的算例中，计算边界处的导数用到一个单边差分格式（公式 ( 15 ) (15) (15)），如下
Γ ? ? ? x ∣ A = D A ? 3 ( 9 ? P ? 8 ? A ? ? E ) (22) \left. \Gamma \frac{\partial \phi}{\partial x} \right |_A = \frac{D^*_A}{3}(9\phi_P-8\phi_A-\phi_E) \tag{22} Γ?x???∣∣∣∣?A?=3DA???(9?P??8?A???E?)(22)
其实就是使用了节点 P 、 E P、E P、E和边界点 A A A的值来近似边界 A A A处的一阶导数，上式可写成
? ? ? x ∣ A = 9 ? P ? 8 ? A ? ? E 3δ x (23) \left. \frac{\partial \phi}{\partial x} \right |_A = \frac{9\phi_P-8\phi_A-\phi_E}{3\ \delta x} \tag{23} ?x???∣∣∣∣?A?=3 δx9?P??8?A???E??(23)
下面用待定系数法来推导一下公式 ( 23 ) (23) (23)所示差分格式的系数及其精度。
先引入待定系数，则
( ? ? ? x ) A = a? A + b? P + c? E (24) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right )_A = a\ \phi_A + b\ \phi_P + c \ \phi_E \end{aligned} \tag{24} (?x???)A?=a ?A?+b ?P?+c ?E??(24)
函数 ? \phi ?在边界 A A A处的泰勒级数展开有
? P = ? A + δ x 2 ( ? ? ? x ) A + 1 2 ( δ x 2 ) 2 ( ? 2 ? ? x 2 ) A + O ( δ x 3 ) ? E = ? A + 3 δ x 2 ( ? ? ? x ) A + 1 2 ( 3 δ x 2 ) 2 ( ? 2 ? ? x 2 ) A + O ( δ x 3 ) (25) \begin{aligned} \phi_P=\phi_A+ \frac{\delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_A+ \frac{1}{2} \left ( \frac{\delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_A +O(\delta x^3) \\\\ \phi_E=\phi_A+ \frac{3\delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_A+ \frac{1}{2} \left ( \frac{3\delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_A +O(\delta x^3) \end{aligned} \tag{25} ?P?=?A?+2δx?(?x???)A?+21?(2δx?)2(?x2?2??)A?+O(δx3)?E?=?A?+23δx?(?x???)A?+21?(23δx?)2(?x2?2??)A?+O(δx3)?(25)
把公式 ( 25 ) (25) (25)代入公式 ( 24 ) (24) (24)中，得
( ? ? ? x ) A = a? A + b [ ? A + δ x 2 ( ? ? ? x ) A + 1 2 ( δ x 2 ) 2 ( ? 2 ? ? x 2 ) A + O ( δ x 3 ) ] + c [ ? A + 3 δ x 2 ( ? ? ? x ) A + 1 2 ( 3 δ x 2 ) 2 ( ? 2 ? ? x 2 ) A + O ( δ x 3 ) ] (26) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right )_A &= a\ \phi_A\\\\ &+b \left[\phi_A+ \frac{\delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_A+ \frac{1}{2} \left ( \frac{\delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_A+O(\delta x^3) \right] \\\\ &+c \left[\phi_A+ \frac{3\delta x}{2} \left( \frac{\partial \phi}{\partial x} \right)_A+ \frac{1}{2} \left ( \frac{3\delta x}{2} \right)^2 \left( \frac{\partial^2\phi}{\partial x^2} \right)_A+O(\delta x^3) \right] \end{aligned} \tag{26} (?x???)A??=a ?A?+b[?A?+2δx?(?x???)A?+21?(2δx?)2(?x2?2??)A?+O(δx3)]+c[?A?+23δx?(?x???)A?+21?(23δx?)2(?x2?2??)A?+O(δx3)]?(26)
然后两边对应系数相等，有
{ a + b + c = 0 δ x 2b + 3 δ x 2c = 1 δ x 2 8 b + 9 δ x 2 8 c = 0 (27) \left \{ \begin{aligned} a + b + c = 0\\\\ \frac{\delta x}{2}\ b + \frac{3\delta x}{2}\ c =1\\\\ \frac{\delta x^2}{8} b + \frac{9\delta x^2}{8} c = 0 \end{aligned} \right. \tag{27} ??????????????????????a+b+c=02δx? b+23δx? c=18δx2?b+89δx2?c=0?(27)
解方程组得
{ a = ? 8 3 δ x b = 9 3 δ x c = ? 1 3 δ x (28) \left \{ \begin{aligned} a&=-\frac{8}{3\delta x}\\\\ b&=\frac{9}{3\delta x}\\\\ c&=-\frac{1}{3\delta x} \end{aligned} \right. \tag{28} ????????????????????????abc?=?3δx8?=3δx9?=?3δx1??(28)
即
( ? ? ? x ) A = ? 8 3 δ x ? A + 9 3 δ x ? P ? 1 3 δ x? E + O ( δ x 2 ) = 9 ? P ? 8 ? A ? ? E 3δ x + O ( δ x 2 ) (29) \begin{aligned} \left( \frac{\partial \phi}{\partial x} \right )_A &=-\frac{8}{3\delta x} \phi_A +\frac{9}{3\delta x} \phi_P -\frac{1}{3\delta x} \ \phi_E +O(\delta x^2)\\\\ &= \frac{9\phi_P-8\phi_A-\phi_E}{3\ \delta x} +O(\delta x^2) \end{aligned} \tag{29} (?x???)A??=?3δx8??A?+3δx9??P??3δx1? ?E?+O(δx2)=3 δx9?P??8?A???E??+O(δx2)?(29)
上式中的截断项 O ( δ x 2 ) O(\delta x^2) O(δx2)是二次方而不是三次方，因为系数 a 、 b 、 c a、b、c a、b、c的分母中都带有 δ x \delta x δx，把系数 a 、 b 、 c a、b、c a、b、c代入公式 ( 26 ) (26) (26)则原来的 O ( δ x 3 ) O(\delta x^3) O(δx3)就变成了 O ( δ x 2 ) O(\delta x^2) O(δx2)。可见，该差分格式是二阶精度的。
参考资料 Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
任玉新. 计算流体力学讲义，2003.