CN115238611A - 一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法 - Google Patents

Abstract

本发明公开了一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，步骤如下：建立对应多相流问题的物理模型，并对整个计算区域进行网格划分；根据问题类型输入计算所需的初始条件以及边界条件；通过相场控制的Cahn–Hilliard方程、以及流场控制的格子玻尔兹曼通量求解器，在有限体积法的基础上进行空间离散，设定全局时间步长，在流动求解过程中使用双时间步推进方法、并利用二阶TVD龙格库塔时间格式进行内迭代伪时间步上的推进计算；在求解过程中不断计算局部时间步并将其应用到双时间步推进方法的迭代过程，直到内迭代残差符合收敛要求为止，得到下一时刻下整个计算域的流场相场数值。本发明能够有效减少多相流问题计算时间，提高计算稳定性。

Description

一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化 方法

技术领域

本发明涉及多相流体计算领域，具体地说是一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法。

背景技术

多相流问题一直是计算流体力学及其他物理学科极为关注的问题，在航空航天、化工生产、能源开发等领域广泛存在，如飞机结冰、石油开采等方面都有其身影，所以，对于多相流问题的研究，始终是科学界关注的焦点，其相关的理论算法具有广阔的发展前景。

为了研究多相流问题，有效模拟诸如不可压缩流体的多相流动，人们提出了很多的数值方法，诸如流体体积法(VOF)、标记单元法(MAC)、扩散界面法(DI)、水平集方法等等。其中，近年来，扩散界面法在处理复杂多相流问题时，具有更强的适应性，得到了显著的普及。扩散界面法作为一种具有较强适应性的方法，其在国内外被广为改进应用。对于现有的扩散界面法主要分为两类。第一类为基于宏观守恒定律的连续介质法，第二类则为基于介观格子玻尔兹曼方程的动理学方法。连续介质方法主要思路为直接求解相场的Cahn–Hilliard方程和流场的Navier-Stokes方程，能够较为容易地应用高阶迎风等格式获得比较精确的模拟结果。而动理学方法基于格子玻尔兹曼方程，由于其考虑了微观粒子的碰撞运动，使得动理学方法能够更好从微观角度出发对流体运动进行模拟。

近年来，在这两种方法的基础上，又诞生了一种多相格子玻尔兹曼通量求解方法(MLBFS)，该方法能够将目前主流两类扩散界面法的优势结合起来，比较精确计算多相流的流场情况。然而，由于许多多相流问题密度比较大、粘性差高，使得在对多相流问题计算时存在稳定性低、计算时间长等问题。对于多相流的计算，需要一种新的优化算法以提高计算效率及稳定性。

发明内容

本发明的目的是针对多相格子玻尔兹曼通量算法目前所存在的问题，提供一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法；该效率优化方法能够减少多相流问题计算时间的优点，适用于大密度比、高粘性比等多相流问题。

本发明的目的是通过以下技术方案解决的：

一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：该效率优化方法的步骤如下：

S1：建立对应多相流问题的物理模型，并对整个计算区域进行网格划分；

S2：根据问题类型输入计算所需的初始条件以及边界条件；

S3：通过相场控制的Cahn–Hilliard方程、以及流场控制的格子玻尔兹曼通量求解器，在有限体积法的基础上进行空间离散，设定全局时间步长，在流动求解过程中使用双时间步推进方法、并利用二阶TVD龙格库塔时间格式进行内迭代伪时间步上的推进计算；

S4：在步骤S3的求解过程中不断计算局部时间步长并将其应用到双时间步推进方法的迭代过程，直到内迭代残差符合预设收敛要求为止，得到下一时刻下整个计算域的流场相场数值。

所述步骤S4中，为了加速内迭代残差收敛，采用隐式残值光顺方法提升局部时间步长以加快收敛速度。

所述的隐式残值光顺方法的公式为：

式(7)中：ε^I、ε^J分别为计算域I、计算域J方向上的松弛系数；

为第一次I方向光顺后(I-1,J)网格中心记录残值；

为第一次I方向光顺后(I,J)网格中心记录残值；

为第一次I方向光顺后(I+1,J)网格中心记录残值；

为未光顺时(I,J)网格中心所记录的初始残值；

为第二次J方向光顺后(I-1,J)网格中心记录最终残值；

为第二次J方向光顺后(I,J)网格中心记录最终残值；

为第二次J方向光顺后(I+1,J)网格中心记录最终残值；

为第一次I方向光顺后(I,J)网格中心记录残值。

所述步骤S1中的物理模型为两相流动模型。

所述步骤S2中的初始条件中的相场的初始分布由相场距离符号函数获得，所选取的相场距离符号函数为：

式(1)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；d则为该控制体到多相流体交界面处的垂直距离，其符号由初始控制体密度决定，对于两相流动，高密度流体处控制体距离为正、低密度流体控制体距离为负；ξ为定义的界面厚度。

所述步骤S2中的边界条件根据研究对象正确选用边界条件。

所述步骤S3中的相场控制的Cahn–Hilliard方程为：

式(2)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；t为时间；

为散度算子；u为流场速度矢量；Γ为质量流率；μ_C代表化学势，且μ_C为：

式(3)中：σ为表面张力系数；ξ为界面厚度。

所述步骤S3中的流场控制的格子玻尔兹曼通量求解器为：

式(4)中：W为宏观守恒量；t为时间；

为散度算子；F为流动通量；S为源项；

式(5)中：p为所在网格格心宏观压强；ρ为所在网格格心宏观压强；u为宏观水平速度；v为宏观竖直速度；c_s为声速，由格子速度模型给出；T为该矩阵的转置；α为所代表的粒子运动方向；e_α代表对应方向上的粒子速度，由格子速度模型给出；

为碰撞后对应方向粒子分布函数；

为碰撞前对应方向粒子分布函数；

为综合分布函数；

为梯度算子；F_s为表面张力；r为对应方向上位移；t为时间；w_α为对应方向上的权重系数，由格子速度模型给出；η＝min(λ,0.025)，λ＝tanh(|C_L-C_R|/(2(C_L+C_R+0.2)))，C_L为对应网格某边左侧控制体的序参数、C_R为对应网格某边右侧控制体的序参数；δ_t为给定粒子运动时间；τ为松弛时间。

非定常流动问题时间推进原始公式为：

式中，Ω为控制体体积，R为残值。

所述步骤S3和步骤S4中的双时间步推进方法的公式为：

式(6)中：Ω为控制体体积；

为伪时间步下的宏观守恒量；t^*为伪时间；R为残值；

为第n步真实时间步下的宏观守恒量；

为第n-1步真实时间步下的宏观守恒量；Δt为全局时间步；R^*为隐式残值光顺后网格中心所记录的最终残值。

本发明相比现有技术有如下优点：

本发明所提供的多相流模拟效率优化方法通过求解相场控制的Cahn–Hilliard方程以及利用格子玻尔兹曼通量求解器计算流场宏观物理量；求解过程中基于双时间步推进方程，在内迭代过程中引入局部时间步长以及隐式残值光顺技术，对多相格子玻尔兹曼通量方法(MLBFS)进行了改进，既保留了原方法在求解多相流问题时所具备的良好的介观特性、又使得在计算大密度比和高粘性比等多相流问题时，能够减少计算时间和提高计算稳定性，增强了MLBFS在多相流体计算领域的适应性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。

附图1为本发明的一个实施例提供的一种基于多相格子玻尔兹曼通量方法(MLBFS)的效率优化方法的方法流程图；

附图2为本发明的一个实施例提供的一种D2Q9格子速度模型的离散速度示意图；

附图3为本发明的一个实施例提供的一种单气泡上升问题的几何模型示意图；

附图4为本发明的一个实施例提供的一种基于低参考速度下的单气泡上升问题的气泡界面位置与原始MLBFS结果对比图；

附图5为本发明的一个实施例提供的一种基于低参考速度下的单气泡上升问题的内迭代步数结果图；

附图6为本发明的一个实施例提供的一种基于高参考速度下的单气泡上升问题的气泡界面位置与原始MLBFS结果对比图；

附图7为本发明的一个实施例提供的一种基于高参考速度下的单气泡上升问题的内迭代步数结果图。

具体实施方式

为了使本技术领域的人员更好地理解本发明中的技术方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

虽然本发明提供了如下述实施例或附图所示的方法操作步骤或装置结构，但基于常规或者无需创造性的劳动在所述方法或装置中可以包括更多或者更少的操作步骤或模块单元。在逻辑性上不存在必要因果关系的步骤或结构中，这些步骤的执行顺序或装置的模块结构不限于本发明实施例或附图所示的执行顺序或模块结构。所述的方法或模块结构的在实际中的装置或终端产品应用时，可以按照实施例或者附图所示的方法或模块结构进行顺序执行或者并行执行(例如并行处理器或者多线程处理的环境、甚至包括分布式处理的实施环境)。

图1为本发明的一个实施例提供的一种基于多相格子玻尔兹曼通量方法(MLBFS)的效率优化方法的方法流程图。本发明的思路是基于双时间步推进方程，在内迭代过程中引入局部时间步长以及隐式残值光顺技术，提高MLBFS对多相流问题的计算效率。

一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法的流程如图1所示，包括：

S1：建立对应多相流问题的物理模型，并对整个计算区域进行网格划分。

S2：根据问题类型输入计算所需的初始条件以及边界条件；初始条件中的相场的初始分布由相场距离符号函数获得，所选取的相场距离符号函数为：

式(1)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；d则为该控制体到多相流体交界面处的垂直距离，其符号由初始控制体密度决定，对于两相流动，高密度流体处控制体距离为正、低密度流体控制体距离为负；ξ为定义的界面厚度。

S3：通过相场控制的Cahn–Hilliard方程、以及流场控制的格子玻尔兹曼通量求解器，在有限体积法的基础上进行空间离散，设定全局时间步长，在流动求解过程中使用双时间步推进方法、并利用二阶TVD龙格库塔时间格式进行内迭代伪时间步上的推进计算；

在利用相场控制的Cahn–Hilliard方程求解相场变化时，其方程为：

式(2)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；t为时间；

为散度算子；u为流场速度矢量；Γ为质量流率；μ_C代表化学势，且μ_C为：

式(3)中：σ为表面张力系数；ξ为界面厚度。

在利用格子玻尔兹曼通量求解器时，其方程为：

式(4)中：W为宏观守恒量；t为时间；

为散度算子；F为流动通量；S为源项；式(4)中的各项的展开式为：

式(5)中：p为所在网格格心宏观压强；ρ为所在网格格心宏观压强；u为宏观水平速度；v为宏观竖直速度；c_s为声速，由格子速度模型给出；T为该矩阵的转置；α为所代表的粒子运动方向；e_α代表对应方向上的粒子速度，由格子速度模型给出；

为碰撞后对应方向粒子分布函数；

为碰撞前对应方向粒子分布函数；

为综合分布函数；

为梯度算子；F_s为表面张力；r为对应方向上位移；t为时间；w_α为对应方向上的权重系数，由格子速度模型给出；η＝min(λ,0.025)，λ＝tanh(|C_L-C_R|/(2(C_L+C_R+0.2)))，C_L为对应网格某边左侧控制体的序参数、C_R为对应网格某边右侧控制体的序参数；δ_t为给定粒子运动时间；τ为松弛时间；由于本发明主要使用如图2所示的D2Q9格子速度模型：故e_α代表对应方向上的粒子速度，

F_s为表面张力，对应系数w₀＝4/9、w₁＝w₂＝w₃＝w₄＝1/9、w₅＝w₆＝w₇＝w₈＝1/36。

计算过程中采用双时间步推进方法的公式以利用内迭代计算优势，为加快计算速度，采用二阶TVD龙格库塔时间格式。双时间步推进方法的公式为：

式(6)中：Ω为控制体体积；

为伪时间步下的宏观守恒量；t^*为伪时间；R为残值；

为第n步真实时间步下的宏观守恒量；

为第n-1步真实时间步下的宏观守恒量；Δt为全局时间步；R^*为隐式残值光顺后网格中心所记录的最终残值。

S4：在步骤S3的求解过程中不断计算局部时间步长并将其应用到双时间步推进方法的迭代过程，直到内迭代残差符合预设收敛要求为止，得到下一时刻下整个计算域的流场相场数值；且为了加速内迭代残差收敛，采用隐式残值光顺方法提升局部时间步长以加快收敛速度。隐式残值光顺方法的公式为：

式(7)中：ε^I、ε^J分别为计算域I、计算域J方向上的松弛系数；

为第一次I方向光顺后(I-1,J)网格中心记录残值；

为第一次I方向光顺后(I,J)网格中心记录残值；

为第一次I方向光顺后(I+1,J)网格中心记录残值；

为未光顺时(I,J)网格中心所记录的初始残值；

为第二次J方向光顺后(I-1,J)网格中心记录最终残值；

为第二次J方向光顺后(I,J)网格中心记录最终残值；

为第二次J方向光顺后(I+1,J)网格中心记录最终残值；

为第一次I方向光顺后(I,J)网格中心记录残值。在本发明给出的实施例中：ε^I＝ε^J＝0.75。

实施例

本发明通过一个具体的实施例来进一步阐述本发明提供的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法及其优越性。

该实施例对大密度比下的单气泡上升问题进行了模拟，几何模型情况说明如图3所示，计算区域网格划分为241×481，参考长度为气泡直径D＝120，整个无量纲计算区域为[0,2D]×[0,4D]，初始气泡被放置在(D,D)位置，上、下两壁面采用无滑移边界条件，而左、右两壁面则采用周期边界条件。

该单气泡上升问题由两个无量纲数决定其雷诺数与厄特沃什数定义方式为：

式(8)中：特征速度

g为重力加速度；σ为表面张力系数；较重流体密度ρ_H＝1，较轻流体密度ρ_L＝0.001；U＝0.0012、Re＝35、Eo＝125，界面厚度ξ＝4，较重流体与较轻流体粘性比μ_H/μ_L＝100，参考时间T＝D/U。

对应的，为保证计算不会发散，采用原始MLBFS进行计算时使用Δt＝1的全局时间步长，而为了说明本发明对可选择的全局时间的影响，使用本发明效率优化算法采用Δt＝10的全局时间步长，计算时间范围为t＝0-7。

对于本实施例中，松弛时间τ可取值为：

式(9)中：δ_t＝δ_x＝0.498，

c＝δ_x/δ_t。

对应的，采用原始MLBFS方法与采用了效率优化后的本发明方法所计算的参考速度为U＝0.0012条件下，t＝5时气泡界面的对比情况图如图4所示。其中，present代表了采用效率优化方法后气泡的界面位置、而original则为原始MLBFS方法所计算的气泡界面位置结果。

所示的效率优化方法对于多相流问题的界面计算结果与原始MLBFS方法对于本实施例的界面计算结果基本保持一致，表明本发明所提出的效率优化方法并不会影响多相流问题的计算准确性，证明了效率优化方法的可行性。

进一步，在利用效率优化算法进行计算时提取每次全局时间步推进时所需的内迭代步数，并将提取后的结果按照时间顺序排列分布，如图5所示。

根据图5所示，内迭代步数的分布较为分散，为了明确本发明对于计算效率优化的有效性，提取了原始MLBFS的内迭代步数与基于低参考速度条件下效率优化算法总步数进行了比对，采用本发明的方法的基于低参考速度下的单气泡上升问题的内迭代总步数与原始MLBFS总步数的对比以及效率提升情况比对结果如表1所示。

算例	方法类型	特征速度U	内迭代总步数	提升效率
					1	原始方法	0.0012	700000	0
2	效率优化方法	0.0012	600516	14.2％

表1基于低参考速度下的单气泡上升问题求解对比

根据表1所示，相比于原始MLBFS计算所需的总步数，效率优化方法计算所需的内迭代总步数，已经从700000步下降到600516步，计算效率提升了14.2％，可见，对于多相流问题的计算，本发明所提出的效率优化算法有效性显著。

在本申请的一个实施例的基础上，在保留其它设置不变的情况下，通过改正特征速度的大小以明确本发明对于多相流计算稳定性的影响以及效率优化的可靠性。

对于修改特征速度，在保留其它初始计算设置不变的条件下，特征速度U＝0.006，利用本发明的效率优化方法进行多相流问题的计算。在该设置条件下，对于原始MLBFS方法，在特征速度超过0.0012时计算容易出不稳定问题从而发散最终计算失败。因此，选择较大的特征速度以验证本发明所提方法在提高稳定性的同时，能够缩小计算周期，从而实现计算效率上的更大突破。

在U＝0.006的设置条件下，使用本发明效率优化算法采用Δt＝10的全局时间步长，计算时间范围为t＝0-7。

对应的，原始MLBFS方法在参考速度为U＝0.0012条件与采用了效率优化后的本发明方法所计算的参考速度为U＝0.006条件下，t＝5时气泡界面的对比情况图如图6所示。其中，present代表了采用效率优化方法后气泡的界面位置，而original则为原始MLBFS方法所计算的气泡界面位置结果。

所示的效率优化方法对于多相流问题的界面计算结果与原始MLBFS方法对于本实施例的界面计算结果基本保持一致，表明本发明所提出的效率优化方法并不会影响多相流问题的计算准确性，证明了效率优化方法的可行性。

进一步，在利用效率优化算法进行计算时提取每次全局时间步推进时所需的内迭代步数，并将提取后的结果按照时间顺序排列分布，如图7所示。

根据图7所示，利用效率优化方法提升特征速度后可以有效计算该类大密度比多相流问题，且内迭代步数的变化幅度相比于图5，有明显下降，可知本发明所提的效率优化方法可以明显提高对于多相流问题计算的稳定性。

为了明确本发明对于计算效率优化的有效性，提取了原始MLBFS的内迭代步数与基于高参考速度条件下的效率优化算法总步数进行了比对，

采用本发明的方法的基于高参考速度下的单气泡上升问题的内迭代总步数与原始MLBFS总步数的对比以及效率提升情况比对结果如表2所示。

算例	方法类型	特征速度U	内迭代总步数	提升效率
					1	原始方法	0.0012	700000	0
2	效率优化方法	0.0060	185845	73.5％

表2基于高参考速度下的单气泡上升问题求解对比

根据表2所示，相比于原始MLBFS计算所需的总步数，效率优化方法计算所需的内迭代总步数，已经从700000步下降到185845步，计算效率提升了73.5％，可见，对于多相流问题的计算，本发明所提出的效率优化算法有效性显著。

结合图1-7与表1-2，相比于原始MLBFS方法，本发明所提出的效率优化方法在双时间法的基础上，通过引入局部时间步长、隐式残值光顺技术，能够大幅度提高全局时间步长，减少计算时间。对于参考速度受到限制的多相流问题，可以较大程度提高参考速度，缩小计算周期，提高计算效率，且计算稳定性也得到了一定程度的提升。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内；本发明未涉及的技术均可通过现有技术加以实现。

Claims (8)

1.一种基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：该效率优化方法的步骤如下：

S1：建立对应多相流问题的物理模型，并对整个计算区域进行网格划分；

S2：根据问题类型输入计算所需的初始条件以及边界条件；

S3：通过相场控制的Cahn–Hilliard方程、以及流场控制的格子玻尔兹曼通量求解器，在有限体积法的基础上进行空间离散，设定全局时间步长，在流动求解过程中使用双时间步推进方法、并利用二阶TVD龙格库塔时间格式进行内迭代伪时间步上的推进计算；

S4：在步骤S3的求解过程中不断计算局部时间步长并将其应用到双时间步推进方法的迭代过程，直到内迭代残差符合预设收敛要求为止，得到下一时刻下整个计算域的流场相场数值。

2.根据权利要求1所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S4中，为了加速内迭代残差收敛，采用隐式残值光顺方法提升局部时间步长以加快收敛速度。

3.根据权利要求2所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述的隐式残值光顺方法的公式为：

式(7)中：ε^I、ε^J分别为计算域I、计算域J方向上的松弛系数；

为第一次I方向光顺后(I-1,J)网格中心记录残值；

为第一次I方向光顺后(I,J)网格中心记录残值；

为第一次I方向光顺后(I+1,J)网格中心记录残值；

为未光顺时(I,J)网格中心所记录的初始残值；

为第二次J方向光顺后(I-1,J)网格中心记录最终残值；

为第二次J方向光顺后(I,J)网格中心记录最终残值；

为第二次J方向光顺后(I+1,J)网格中心记录最终残值；

为第一次I方向光顺后(I,J)网格中心记录残值。

4.根据权利要求1-3任一所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S1中的物理模型为两相流动模型。

5.根据权利要求1-3任一所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S2中的初始条件中的相场的初始分布由相场距离符号函数获得，所选取的相场距离符号函数为：

式(1)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；d则为该控制体到多相流体交界面处的垂直距离，其符号由初始控制体密度决定，对于两相流动，高密度流体处控制体距离为正、低密度流体控制体距离为负；ξ为定义的界面厚度。

6.根据权利要求1-3任一所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S3中的相场控制的Cahn–Hilliard方程为：

式(2)中：C代表相场内的序参数，即对应控制体内较重质量流体的体积分数；t为时间；

为散度算子；u为流场速度矢量；Γ为质量流率；μ_C代表化学势，且μ_C为：

式(3)中：σ为表面张力系数；ξ为界面厚度。

7.根据权利要求1-3任一所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S3中的流场控制的格子玻尔兹曼通量求解器为：

式(4)中：W为宏观守恒量；t为时间；

为散度算子；F为流动通量；S为源项；

式(5)中：p为所在网格格心宏观压强；ρ为所在网格格心宏观压强；u为宏观水平速度；v为宏观竖直速度；c_s为声速，由格子速度模型给出；T为该矩阵的转置；α为所代表的粒子运动方向；e_α代表对应方向上的粒子速度，由格子速度模型给出；

为碰撞后对应方向粒子分布函数；

为碰撞前对应方向粒子分布函数；

为综合分布函数；

为梯度算子；F_s为表面张力；r为对应方向上位移；t为时间；w_α为对应方向上的权重系数，由格子速度模型给出；η＝min(λ,0.025)，λ＝tanh(|C_L-C_R|/(2(C_L+C_R+0.2)))，C_L为对应网格某边左侧控制体的序参数、C_R为对应网格某边右侧控制体的序参数；δ_t为给定粒子运动时间；τ为松弛时间。

8.根据权利要求1-3任一所述的基于多相格子玻尔兹曼通量方法的多相流模拟效率优化方法，其特征在于：所述步骤S3和步骤S4中的双时间步推进方法的公式为：

式(6)中：Ω为控制体体积；

为伪时间步下的宏观守恒量；t^*为伪时间；R为残值；

为第n步真实时间步下的宏观守恒量；

为第n-1步真实时间步下的宏观守恒量；Δt为全局时间步；R^*为隐式残值光顺后网格中心所记录的最终残值。

Images