本发明实施例提供一种两相流动界面演化模拟方法及系统，包括：建立关于该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程，引入稳定项，空间离散后得到该两相流相场模型的半离散系统；采用显式积分因子方法进行时间离散，得到时空全离散格式；将显式积分因子方法中的指数函数进行合理逼近，得到无条件保持该两相流相场模型极值原理与质量守恒的显式单步时空全离散格式；得到下一时间层上的界面的近似值，依次迭代得到计算区域内终止时刻两相流相场界面的模拟结果。基于时间离散方法的改进处理，特别是稳定化参数的选取以及对指数函数的逼近，使得时间离散格式可以达到四阶无条件保持极值原理与质量守恒。
一种两相流动界面演化模拟方法及系统
技术领域
[0001] 本发明涉及两相流动仿真领域，具体涉及一种两相流动界面演化模拟方法及系统。
背景技术
[0002] 两相流动广泛存在于自然界与现代工业生产过程中，与人类的生产与生活密切相关。基于计算流体力学的工业模拟的迅猛发展趋势对复杂两相流动仿真的需求，对两相流模型及算法的计算效果提出了越来越高的要求。针对相场Allen‑Cahn方程及其保持质量守恒律的修正方程，当前大多数格式在保持极值原理时对时间步长有严格的要求，而高阶的积分因子Runge‑Kutta(IFRK)格式在时间步长较大的时候又会出现指数衰减现象，这给快速求解Allen‑Cahn方程、高效捕捉两相界面带来困难。
发明内容
[0003] 本发明实施例提供一种两相流动界面演化模拟方法及系统，基于时间离散方法的改进处理，特别是稳定化参数的选取以及对指数函数的逼近，使得时间离散格式可以达到四阶无条件保持极值原理与质量守恒。
[0004] 为达上述目的，一方面，本发明实施例提供一种两相流动界面演化模拟方法，包括如下步骤：
[0005] 步骤一：建立关于该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程，并在守恒型相场Allen‑Cahn方程中引入稳定项，之后采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统；
[0006] 步骤二：采用显式积分因子Runge‑Kutta方法对该两相流相场模型的半离散系统进行时间离散，得到该两相流的时空全离散格式；
[0007] 步骤三：针对该两相流相场模型的时空全离散格式，将显式积分因子Runge‑Kutta方法中的指数函数进行合理逼近，得到无条件保持该两相流相场模型极值原理与质量守恒的显式单步时空全离散格式；
[0008] 步骤四：根据该两相流相场模型的无条件保持极值原理与质量守恒的显式单步全离散格式，得到该两相流相场模型在下一时间层上的界面的近似值，依次迭代，得到计算区域内终止时刻两相流相场界面的模拟结果，通过终止时刻两相流相场界面的模拟结果表示该两相流界面形态。
[0009] 另一方面，本发明实施例提供一种两相流动界面演化模拟系统，包括：
[0010] 处理器；以及
[0011] 被安排成存储计算机可执行指令的存储器，所述可执行指令在被执行时使所述处理器执行以下操作：
[0012] 步骤一：建立关于该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程，并在守恒型相场Allen‑Cahn方程中引入稳定项，之后采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统；
[0013] 步骤二：采用显式积分因子Runge‑Kutta方法对该两相流相场模型的半离散系统进行时间离散，得到该两相流的时空全离散格式；
[0014] 步骤三：针对该两相流相场模型的时空全离散格式，将显式积分因子Runge‑Kutta方法中的指数函数进行合理逼近，得到无条件保持该两相流相场模型极值原理与质量守恒的显式单步时空全离散格式；
[0015] 步骤四：根据该两相流相场模型的无条件保持极值原理与质量守恒的显式单步全离散格式，得到该两相流相场模型在下一时间层上的界面的近似值，依次迭代，得到计算区域内终止时刻两相流相场界面的模拟结果，通过终止时刻两相流相场界面的模拟结果表示该两相流界面形态。
[0016] 上述技术方案具有如下有益效果：基于时间离散方法的改进处理，特别是稳定化参数的选取以及对指数函数的逼近，使得时间离散格式可以达到四阶无条件保持极值原理与质量守恒。
附图说明
[0017] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0018] 图1是本发明实施例的两相流动界面演化模拟方法的流程图；
[0019] 图2是本发明实施例的基于四阶isIFRK方法的三个凹八面体融合问题模拟结果示意图；
[0020] 图3是本发明实施例的使用sIFRK与isIFRK计算的数值解对极值原理和质量守恒律的保持。
具体实施方式
[0021] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0022] 如图1所示，结合本发明的实施例，提供一种两相流动界面演化模拟方法，包括如下步骤：
[0023] 步骤一：建立关于该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程，并在守恒型相场Allen‑Cahn方程中引入稳定项，之后采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统；
[0024] 步骤二：采用显式积分因子Runge‑Kutta方法对该两相流相场模型的半离散系统进行时间离散，得到该两相流的时空全离散格式；
[0025] 步骤三：针对该两相流相场模型的时空全离散格式，将显式积分因子Runge‑Kutta方法中的指数函数进行合理逼近，得到无条件保持该两相流相场模型极值原理与质量守恒的显式单步时空全离散格式；
[0026] 步骤四：根据该两相流相场模型的无条件保持极值原理与质量守恒的显式单步全离散格式，得到该两相流相场模型在下一时间层上的界面的近似值，依次迭代，得到计算区域内终止时刻两相流相场界面的模拟结果，通过终止时刻两相流相场界面的模拟结果表示该两相流界面形态。
[0027] 优选地，所述步骤一，具体包括：
[0028] 将该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程表示为：
[0029]
[0030] 其中，x表示空间方向坐标，Ω表示一维、二维或三维区域，t表示时间坐标，T表示终止时刻，u表示相场变量，ut表示相场变量关于时间的导数，Δ表示二阶拉普拉斯算子，是带有拉格朗日乘子的非线性项， 表示为 f(u)为势函数的负导数，即f(u)＝‑F'(u)；u0(x)表示u在t＝0时刻的解，u0(x)表示相场的初始值；相场变量u(x)＝1表示在x处为相一，u(x)＝‑1表示在x处为相二，‑1＜u(x)＜1表示x处相一、相二均存在，该处为两相的接触面。
[0031] 当势函数采用Ginzburg‑Landau势函数F(u)： 时，f(u)表示为：
[0032] f(u)＝u‑u3；
[0033] 对于拉格朗日乘子，将其选择为第一拉格朗日乘子RSLM或者第二拉格朗日乘子BBLM，其中：
[0034] 第一拉格朗日乘子RSLM表示为：
[0035] 第二拉格朗日函数BBLM表示为：
[0036] 在式(3)中，λ0(t)表示是式(3)函数的记号，无物理意义；
[0037] 将方程(1)内的拉格朗日乘子选择为RSLM或者BBLM时，该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程(1)满足质量守恒律且该两相流在周期边界条件下界面演化的守恒型相场
Allen‑Cahn方程(1)具有极值原理；
[0038] 将稳定项表示为κu；在两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程(1)中引入稳定项κu得到其等价方程：
[0039]
[0040] 其中，参数κ需要满足如下条件：
[0041]
[0042] 其中，f'(ξ)表示函数f(ξ)关于ξ的一阶导数，ξ表示绝对值小于β的任意一个变量；
[0043] 其中，所述采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统，具体包括：
[0044] 针对两相流相场在一维空间状态下演化时，取Ω＝[xl,xr]，xl为区间左端点，xr为区间右端点，进行等距剖分，采用空间步长h，空间网格节点数为 令u＝[u0,T T
u1,…,uN‑1] 为半离散后的解，其中，表示向量的转置，二阶拉普拉斯算子Δ离散后的二阶中心微分矩阵为：
[0045]N×N 2 2
[0046] 令I∈R 为单位矩阵，记L＝∈ D，Lκ＝∈D‑κI， 则得到在一维空间状态下该两相流相场模型的半离散系统：
[0047] u＝Lκu+Nκ(u) (5)[0048] 在式(5)中，Lκ的下标表κ示稳定项内的参数κ；
[0049] 当κ满足条件(4)时，通过对(5)进行计算得到 且在一维空间状态下该两相流相场模型的半离散系统(5)满足质量守恒律：
[0050]
[0051] 其中，1＝[1,…,1]T∈RN，
[0052] 优选地，所述步骤二，具体包括：
[0053] 对在一维空间状态下该两相流相场模型的半离散系统(5)应用显式积分因子Runge‑Kutta方法，得到在一维空间状态下该两相流相场模型的时空全离散格式sIFRK：
[0054]
[0055] 其中，n表示第n层时间步，ci表示Runge‑Kutta方法的节点，ai,j表示Runge‑Kutta方法的系数，τ表示时间步长，s表示Runge‑Kutta方法的级数。
[0056] 优选地，所述步骤三，具体包括：
[0057] 对在一维空间状态下该两相流相场模型的时空全离散格式sIFRK(6)中的指数函数 采用递归方式进行逼近，令：
[0058]
[0059] 根据方程(7)对该两相流相场模型的时空全离散格式sIFRK(6)进行改进，得到该两相流相场模型改进的显式单步时空全离散格式isIFRK：
[0060]
[0061] 当isIFRK(8)的系数满足0≤cj≤ci≤1，ai,j≥0， 以及当参数κ的选取条件满足方程(4)时，对于任意的时间步长τ＞0时，isIFRK(8)均无条件保持极值原理和离散的质量守恒律
[0062] 优选地，所述步骤四，具体包括：
[0063] 根据该两相流相场模型的显式单步时空全离散格式isIFRK(8)具有关于时间方向迭代的特性，以及根据isIFRK(8)具有极值原理的特征；得到：
[0064] 当isIFRK(8)的初始值满足 时，通过迭代isIFRK(8)求出的下一时间层n+1 0
上的近似值，所述下一时间层上的近似值满足 和＝；
[0065] 依次迭代计算得到在计算区域内终止时刻两相流相场的模拟结果，将终止时刻两相流相场的模拟结果作为该两相流界面形态。
[0066] 优选地，在步骤一中，
[0067] 所述该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程(1)具有极值原理是指:
[0068] 存在常数β，若相场的初始值u0(x)满足 则两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程(1)在之后连续的每一时刻的解均满足当采用Ginzburg‑Landau势函数时，若选择第一拉格朗日乘子RSLM，则
若选择第二拉格朗日乘子BBLM，则β＝1，此时若相场变量u满足极值原理，则表示相场具有物理意义；
[0069] 其中， 表示最大模范数。
[0070] 优选地，在步骤三中，当isIFRK(8)的系数满足0≤cj≤ci≤1，ai,j≥0，以及当参数κ的选取条件满足方程(4)时，对于任意的时间步长τ＞0时,isIFRK(8)均无条件保持极值原理；针对isIFRK(8)对于任意的时间步长τ＞0均无条件保持极值原理证明如下：
[0071] 根据矩阵L为对角元非正的对角占优矩阵，得知 对格式isIFRK(8)两边同时取无穷范数 得到：
[0072]
[0073] 其中，i＝1,…,s，因此得到
[0074] 优选地，在步骤三中，
[0075] 当isIFRK(8)的系数满足0≤cj≤ci≤1，ai,j≥0， 以及当参数κ的选取条件满足方程(4)时，对于任意的时间步长τ＞0时，isIFRK(8)均无条件保持离散的质量守恒律 针对isIFRK(8)对于任意的时间步长τ＞0均无条件保持质量守
恒律证明如下：
[0076] 根据λ＝0为矩阵L对应特征向量1的特征值，得知＝；以及，根据＝κ，对isIFRK(8)两边同时与向量1作内积，得到：
[0077]
[0078] 其中，i＝1,…,s；因此得到＝≡const。
[0079] 得到满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法包括：
[0080]
[0081] 其中，s表示Runge‑Kutta方法的级数，p表示阶数。
[0082] 优选地，所述步骤一，所述采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统，具体包括：
[0083] 对于带有周期边界的二维空间状态和三维空间状态，记x方向的空间网格节点数为Nx，y方向的空间网格节点数为Ny，z方向的空间网格节点数为Nz，单位矩阵Ix：单位矩阵 单位矩阵 单位矩阵 单位矩阵
[0084] 通过使用二阶中心差分方法离散二阶拉普拉斯算子Δ得到二维情况下的微分矩阵： 以及通过使用二阶中心差分方法离散二阶拉普拉斯算子Δ得到三维情况下的微分矩阵：
[0085] 其中， 为克罗内克积，记L2,κ＝∈2D2‑κI2，L3,κ＝∈2D3‑κI3时，得到二维空间状态下的半离散系统：u＝L2,κu+Nκ(u)，以及得到三维空间状态下的半离散系统：u＝L3,κu+Nκ(u)。
[0086] 结合本发明的实施例，还提供一种两相流动界面演化模拟系统，包括：
[0087] 处理器；以及
[0088] 被安排成存储计算机可执行指令的存储器，所述可执行指令在被执行时使所述处理器执行以下操作：
[0089] 步骤一：建立关于该两相流在周期边界条件下界面演化的守恒型相场Allen‑Cahn方程，并在守恒型相场Allen‑Cahn方程中引入稳定项，之后采用二阶中心差分格式对引入稳定项的守恒型相场Allen‑Cahn方程进行空间离散，得到该两相流相场模型的半离散系统；
[0090] 步骤二：采用显式积分因子Runge‑Kutta方法对该两相流相场模型的半离散系统进行时间离散，得到该两相流的时空全离散格式；
[0091] 步骤三：针对该两相流相场模型的时空全离散格式，将显式积分因子Runge‑Kutta方法中的指数函数进行合理逼近，得到无条件保持该两相流相场模型极值原理与质量守恒的显式单步时空全离散格式；
[0092] 步骤四：根据该两相流相场模型的无条件保持极值原理与质量守恒的显式单步全离散格式，得到该两相流相场模型在下一时间层上的界面的近似值，依次迭代，得到计算区域内终止时刻两相流相场界面的模拟结果，通过终止时刻两相流相场界面的模拟结果表示该两相流界面形态。
[0093] 下面结合具体的应用实例对本发明实施例上述技术方案进行详细说明，实施过程中没有介绍到的技术细节，可以参考前文的相关描述。
[0094] 本发明为一种两相流动界面演化模拟方法，也就是一种高阶无条件保持守恒型相场Allen‑Cahn方程极值原理与质量守恒的显式单步法的构造，可以提高守恒型Allen‑Cahn方程模拟效率和精度的计算方法。通过在带有周期边界条件的守恒型相场Allen‑Cahn方程中引入稳定项，然后使用中心差分格式进行空间离散，然后给出了一类时间离散方法的改进处理，特别是稳定化参数的选取以及对指数函数的逼近，使得时间离散格式可以达到四阶无条件保持极值原理与质量守恒。本发明的优点在于时间步长可以任意选择并保持系统的极值原理与质量守恒。通过实例展示了这个新算法在数值模拟中的高效、高精度及保极值与质量守恒等优势。本发明还可以与计算流体力学问题相结合，适用于结合相场Allen‑Cahn方程的数值模拟，最终能够推广到基于相场方法模拟的大规模科学计算问题。
[0095] 1.一种高阶无条件保持相场Allen‑Cahn方程极值原理与质量守恒的显式单步法，包括如下步骤：
[0096] 步骤一：在周期边界边界条件下的守恒型相场Allen‑Cahn方程中引入稳定项，采用二阶中心差分格式进行空间离散，得到半离散系统；
[0097] 步骤二：针对半离散系统采用显式积分因子Runge‑Kutta方法进行时间离散，得到时空全离散系统；
[0098] 步骤三：针对积分因子Runge‑Kutta方法中的指数函数进行合理逼近，得到无条件保持极值原理与质量守恒的显式单步法；
[0099] 步骤四：根据时空全离散无条件保持极值原理与质量守恒的显式单步法，得到下一时间层上的近似值，依次迭代，得到计算区域内终止时刻相场的模拟结果。
[0100] 2.所述步骤一中，守恒型相场Allen‑Cahn方程为
[0101]
[0102] 其中，Ω可以为一维、二维或三维区域, 是带有拉格朗日乘子的非线性项，f(u)＝‑F'(u)为势函数的负导数。当采用Ginzburg‑Landau势函数时有：
[0103]
[0104] 当拉格朗日乘子选择为：
[0105] 第一拉格朗日乘子RSLM表示为：
[0106] 或
[0107] 第二拉格朗日函数BBLM表示为：
[0108] 时方程(1)满足质量守恒律 且具有极值原理：即存在常数β,若初始值满足 则以后每一时刻的解都满足
其中，每一时刻是在连续情况下是任一时刻，也是在离散情况下是时间层。
[0109] 在守恒型相场Allen‑Cahn方程中引入稳定项κu可以得到等价方程：
[0110]
[0111] 其中参数κ的需要满足
[0112]
[0113] 从式(4)可以看出，使用不同拉格朗日乘子时，参数κ(即：kappa)的选择也不同。
[0114] 以一维问题为例，取Ω＝[xl,xr]，进行等距剖分，采用空间步长h，空间网格节点数T为 令u＝[u0,u1,…,uN‑1]为半离散后的解，拉普拉斯算子Δ离散后的二阶中心微分矩阵为
[0115]2 2
[0116] 令 为单位矩阵，记L＝∈D，Lκ＝∈D‑κI， 则空间半离散格式为：
[0117] u＝Lκu+Nκ(u) (5)[0118] 可以证明 且半离散方程(5)具有质量守恒律：
[0119]
[0120] 其中 二维和三维问题的离散方法可以类似给出。
[0121] 3.所述步骤二中，对稳定化的半离散系统(5)应用显式积分因子Runge‑Kutta(integrating factor Runge‑Kutta,IFRK)方法得到：
[0122]
[0123] 记格式(6)为stabilized IFRK(sIFRK)；
[0124] 其中，稳定化指(5)中带有稳定项。
[0125] 4.所述步骤三中，针对格式(6)中的指数函数 采用递归方式进行逼近，令[0126]
[0127] 然后得到改进的积分因子Runge‑Kutta方法(improved sIFRK,isIFRK)：
[0128]
[0129] 当Runge‑Kutta方法的系数满足0≤cj≤ci≤1，ai,j≥0， 可以证明当参数κ的选取满足步骤一中的条件(4)时，格式(8)对于任意的时间步长τ＞0都可以无条件保持极值原理和离散的质量守恒律 无条件即是指不依赖于时间步长，对于任意的时间步长都满足。极值原理证明如下：
τL
[0130] 由矩阵L为对角元非正的对角占优矩阵，知‖e ‖∞≤1， 对格式(8)两边同时取无穷范数 得
[0131]
[0132] 其中，i＝1,…,s；因此得到
[0133] 质量守恒证明如下：
[0134] 由λ＝0为矩阵L对应特征向量1的特征值，知＝，又＝κ，对格式(8)两边同时与向量1作内积，可得
[0135]
[0136] 其中，i＝1,…,s；因此得到＝≡const。
[0137] 满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法包括但不限于：
[0138]
[0139] 表1.满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法
[0140] 5.所述步骤四中，时空全离散无条件保持极值原理与质量守恒的显式单步法为关于时间方向的迭代公式，若初始值满足 通过迭代公式(8)求出的下一时间层的n+1 0
值满足 和＝，依次计算得到终止时刻计算区域内的相场的模拟
结果。
[0141] 为优化上述技术方案，采取的具体措施还包括：
[0142] 所述步骤一中，守恒型相场Allen‑Cahn方程为：
[0143]
[0144] 其中x表示空间方向坐标，Ω表示一维区域、二维区域或三维区域，t表示时间坐标，T表示终止时刻，u表示相场变量，ut表示相场变量关于时间的导数，Δ表示二阶拉普拉斯算子， 是带有拉格朗日乘子的非线性项， 表示为 f(u)＝‑F'(u)，f(u)为势函数的负导数；u(x,0)表示u在初始时刻的解，u0(x)表示相场的初始值；
[0145] 当针对F(u)采用Ginzburg‑Landau势函数时有：
[0146]
[0147] 将拉格朗日乘子选择为：
[0148] 第一拉格朗日乘子RSLM表示为：
[0149] 或
[0150] 第二拉格朗日函数BBLM表示为：
[0151] 当将拉格朗日乘子选择为(2)和(3)时，方程(1)满足质量守恒律且具有极值原理：即存在常数β,若初始值满足
则以后每一时刻的解都满足
[0152] 在守恒型相场Allen‑Cahn方程中引入稳定项κu可以得到等价方程：
[0153]
[0154] 其中，参数κ的需要满足：
[0155]
[0156] 以一维空间状态问题为例，取Ω＝[xl,xr]，进行等距剖分，采用空间步长h，空间网T格节点数为 令u＝[u0,u1,…,uN‑1]为半离散后的解，拉普拉斯算子Δ离散后的二阶中心微分矩阵为：
[0157]N×N 2 2
[0158] 令I∈R 为单位矩阵，记L＝∈ D，Lκ＝∈D‑κI， 则空间半离散格式为：
[0159] u＝Lκu+Nκ(u) (5)[0160] 可以证明 且半离散方程(5)具有质量守恒律：
[0161]
[0162] 其中 二维和三维问题的离散方法可以类似给出。
[0163] 所述步骤二中，对稳定化的空间半离散系统(5)应用显式积分因子Runge‑Kutta(integrating factor Runge‑Kutta,IFRK)方法得到时空全离散系统，所述时空全离散系统表示为:
[0164]
[0165] 记格式(6)为stabilized IFRK(sIFRK)。
[0166] 所述步骤三中，针对格式(6)中的指数函数 采用递归方式进行逼近，令[0167]
[0168] 然后得到改进的积分因子Runge‑Kutta方法(improved sIFRK,isIFRK)：
[0169]
[0170] 当Runge‑Kutta方法的系数满足0≤cj≤ci≤1，ai,j≥0， 可以证明当参数κ的选取满足步骤一中的条件(4)时，格式(8)对于任意的时间步长τ＞0都可以保持极值原理和离散的质量守恒律 格式(8)保持极值原理的证明如下：
[0171] 由矩阵L为对角元非正的对角占优矩阵，知‖eτL‖∞≤1， 对格式(8)两边同时取无穷范数 得
[0172]
[0173] 其中，i＝1,…,s；因此得到
[0174] 质量守恒证明如下：
[0175] 由λ＝0为矩阵L对应特征向量1的特征值，知〈eτLu,1>＝〈u,1>，又＝κ，对格式(8)两边同时与向量1作内积，可得
[0176]
[0177] 其中，i＝1,…,s；因此得到＝≡const
[0178] 满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法包括但不限于：
[0179]
[0180] 表1.满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法
[0181] 所述步骤四中，时空全离散无条件保持极值原理与质量守恒的显式单步法为关于时间方向的迭代公式，若初始值满足 通过迭代公式(8)求出的下一时间层的值n+1 0
满足 和＝，依次计算得到终止时刻计算区域内的相场的模拟结
果。
[0182] 本发明的有益效果是：在现有技术中，计算两相流常用的界面捕捉方法包括：MAC(Mark and Cell Method)方法、VOF(Volume of Fluid Method)方法以及相场(Phase Field Method)方法。1995年，Antannovskii和Jacqmin提出相场方法，随后被用于两相流动的数值模拟领域。相场方法使用连续、变化的相场序参数来描述流体界面的运动，因此可以较为直接的模拟不同相之间的复杂变化过程，已被证明是描述多相流体动力学的有效方法之一。Allen‑Cahn模型及Cahn‑Hillard模型是相场方法中比较常用的两种模型。Allen‑Cahn模型不能保证相变量所描述的物理量在整个区域上积分后守恒(如体积守恒、质量守恒)，而Cahn‑Hillard方程却可以严格地保证这一守恒性。但是Cahn‑Hillard方程中空间算子的阶数要高于Allen‑Cahn方程中空间算子的阶数，这导致数值求解Cahn‑Hilard方程需要的时间成本要高于求解Allen‑Cahn方程。
[0183] 相比于现有技术，通过在经典的Allen‑Cahn模型后添加拉格朗日乘子即可以得到具有守恒性质的修正Allen‑Cahn模型，并且这一模型满足解有界和能量耗散性质。本发明的isIFRK格式(8)通过对sIFRK格式中的指数函数进行合理逼近，消去了sIFRK格式带来的指数衰减或增长现象(由格式(8)保持极值原理的证明可知)，可以无条件保持极值原理和质量守恒，在时间步进中可以允许任意大的时间步长，并保证任意时刻的数值解具有物理意义；相比于已有的数值格式，isIFRK格式由于取消了Runge‑Kutta系数中对增长函数的要求，使得无条件保持极值原理的数值格式最高可以达到4阶。因此该isIFRK格式构造更简单，更易于应用到高维问题的数值模拟中去。所以在开发数值格式时，开发可以保持相场方程极值原理的能量稳定格式。同时避免指数衰减现象，使得格式在时间步长较大时也能给出较为精确的数值解，这对提高两相流动问题的界面捕捉效率具有重要意义。
[0184] 另外，本发明提供的基于改进积分因子格式的高阶无条件保持守恒型相场Allen‑Cahn方程极值原理和质量守恒的显式单步法，可以针对一维、二维和三维守恒型Allen‑Cahn方程进行数值模拟，该模拟方法是显式、高阶的，可以大幅度提高计算效率和计算精度，并且对于任意的时间步长都可以保持系统的极值原理与质量守恒。
[0185] 如图1所示，该模拟方法具体包括如下步骤：
[0186] 一、针对原方程添加稳定项，采用二阶中心差分格式建立空间半离散系统。
[0187] 在三维空间，考虑RSLM类型的守恒型Allen‑Cahn方程
[0188]
[0189] 其中参数∈＝0.01,带有拉格朗日乘子的非线性项 f(u)＝u‑u3,极值为 通过添加稳定项κu，κ≥max|ξ|≤β|f'(ξ)|＝3，令x，y，
z三个方向的网格数均为128，应用中心差分格式得到空间半离散系统
[0190] u＝Lκu+Nκ(u) (1)[0191] 二、建立半离散系统的稳定化积分因子Runge‑Kutta方法。
[0192] 对稳定化的半离散系统应用显式积分因子Runge‑Kutta方法得到sIFRK格式[0193]
[0194] 三、对指数函数进行合理逼近，建立改进的积分因子Runge‑Kutta方法。
[0195] 针对格式(2)中的指数函数 采用递归方式进行逼近，令
[0196]
[0197] 然后得到改进的积分因子Runge‑Kutta方法：
[0198]
[0199] 要求格式(3)中的Runge‑Kutta系数满足0≤cj≤ci≤1，ai,j≥0， 几种一至四阶的Runge‑Kutta方法包括但不限于：
[0200]
[0201] 其中，s表示Runge‑Kutta方法的级数，p表示阶数。
[0202] 表1.满足0≤cj≤ci≤1，ai,j≥0， 的Runge‑Kutta方法
[0203] 四、模拟计算与结果输出
[0204] 基于上述改进的单步法编写计算机程序，形成两相流动界面演化模拟系统，可以获得并导出相场及相关物理量的演化。例如，采用第一拉格朗日乘子RSLM时，基于四阶isIFRK方法的三个凹八面体融合问题模拟结果如图2所示；使用sIFRK与isIFRK计算的数值解对极值原理和质量守恒律的保持如图3所示。
[0205] 应该明白，公开的过程中的步骤的特定顺序或层次是示例性方法的实例。基于设计偏好，应该理解，过程中的步骤的特定顺序或层次可以在不脱离本公开的保护范围的情况下得到重新安排。所附的方法权利要求以示例性的顺序给出了各种步骤的要素，并且不是要限于所述的特定顺序或层次。
[0206] 在上述的详细描述中，各种特征一起组合在单个的实施方案中，以简化本公开。不应该将这种公开方法解释为反映了这样的意图，即，所要求保护的主题的实施方案需要比清楚地在每个权利要求中所陈述的特征更多的特征。相反，如所附的权利要求书所反映的那样，本发明处于比所公开的单个实施方案的全部特征少的状态。因此，所附的权利要求书特此清楚地被并入详细描述中，其中每项权利要求独自作为本发明单独的优选实施方案。
[0207] 为使本领域内的任何技术人员能够实现或者使用本发明，上面对所公开实施例进行了描述。对于本领域技术人员来说；这些实施例的各种修改方式都是显而易见的，并且本文定义的一般原理也可以在不脱离本公开的精神和保护范围的基础上适用于其它实施例。因此，本公开并不限于本文给出的实施例，而是与本申请公开的原理和新颖性特征的最广范围相一致。
[0208] 上文的描述包括一个或多个实施例的举例。当然，为了描述上述实施例而描述部件或方法的所有可能的结合是不可能的，但是本领域普通技术人员应该认识到，各个实施例可以做进一步的组合和排列。因此，本文中描述的实施例旨在涵盖落入所附权利要求书的保护范围内的所有这样的改变、修改和变型。此外，就说明书或权利要求书中使用的术语“包含”，该词的涵盖方式类似于术语“包括”，就如同“包括，”在权利要求中用作衔接词所解释的那样。此外，使用在权利要求书的说明书中的任何一个术语“或者”是要表示“非排它性的或者”。
[0209] 本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrative logical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrative components)，单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
[0210] 本发明实施例中所描述的各种说明性的逻辑块，或单元都可以通过通用处理器，数字信号处理器，专用集成电路(ASIC)，现场可编程门阵列或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。
[0211] 本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD‑ROM或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。
[0212] 在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD‑ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电缆、双绞线、数字用户线(DSL)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。
[0213] 以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
