一种多介质流体界面运动的仿真流程设计方法与流程

本发明属于计算流体力学技术领域，涉及一种多介质流体界面运动的仿真流程设计方法。

背景技术：

随着计算能力的提升，高置信度的数值模拟作用日益增强。传统的cfd软件算法上已经日臻成熟，但是对于多介质流体界面附近的数值处理是一个挑战性难题。这里所谓的界面，是指两种或两种以上不相掺混的介质具有任意拓扑结构的分界面，或者指流体的自由面，通常具有复杂变形、强间断、强非线性以及多物理过程强耦合等特点，广泛存在于航空、航天、航海、核能、材料、天体力学和生物医疗等诸多领域。

由于界面两侧不同流体的物性参数差异很大，如果将模拟单介质流体的数值格式直接应用于模拟多介质流体问题，数值不稳定就会出现，界面附近会出现非物理振荡，甚至使得计算难以进行。有些方法，比如采用混合模型描述界面条件，需要在混合区域里额外构造特殊的状态方程或数值扩散抑制技术，实现起来非常复杂；或者将界面网格切割成包含不同流体的子单元，界面附近小网格问题处理起来非常困难，容易引起数值不稳定。如果将界面看成一类特殊的边界来处理，在界面附近的虚拟流场定义边界条件，可以避免界面附近的小网格问题，也容易推广到多维。但是，上述技术通常只能解决这种单一的多介质界面问题，目前没有设计成求解模块应用于一般的cfd软件，通用性不强。如果需要在cfd软件中实现多介质流体界面问题的模拟，往往需要较大的修改量，程序架构会发生改变，也容易影响cfd软件已有的计算功能和计算效率。

针对以上问题，本发明旨在提出一种多介质流体界面运动的仿真流程设计方法，能够应用于传统的cfd软件。在不影响该软件已有功能和计算效率的基础上，通过添加几个模块和极少量的代码修改实现多介质流体界面问题的数值模拟。该方法实现简单，具有重要的应用价值。

技术实现要素：

本发明提出的一种多介质流体界面运动的仿真流程设计方法，实现对多介质流体界面复杂变形问题的数值模拟，具有极强的通用性，且代码修改量小、内存需求低。

本发明的技术解决方案是：

提供一种多介质流体界面运动的仿真流程设计方法，包括：

(1)输入介质总数nmed，多介质流场仿真所需参数，计算网格信息、边界条件类型；

(2)计算网格单元几何信息、并初始化符号距离函数φi,j,k以及初始化流场qi,j,k，其中i,j,k表示网格节点；

(3)计算时间步长δt；

(4)根据输入的边界条件类型对边界条件赋值；令m＝1；

(5)确定介质m在n时刻的计算区域

(6)建立水平集方程，并求解n+1时刻的符号距离函数

(7)根据n+1时刻的符号距离函数得出介质m在n+1时刻的流场区域进而获得网格节点与介质m的对应关系标识

(8)更新介质m在n时刻的界面条件表示介质m的界面外虚拟流体区域；

(9)计算介质m在n+1时刻的流场判断m是否等于nmed，如果等于则进入步骤(10)，否则令m＝m+1，返回步骤(5)；

(10)计算n+1时刻整个计算域的流场；

(11)判断是否到达结果输出时间，如果到达则进入步骤(12)，否则返回步骤(3)；

(12)输出n+1时刻流场计算结果

(13)输出n+1时刻符号距离函数计算结果

(14)判断时间推进是否结束，如果结束则完成多介质流体界面运动仿真；否则返回(3)。

优选的，符号距离函数表示每个网格节点到界面的符号距离，符号为正值的节点代表一种流体，符号为负值的节点代表另一种流体，符号为零的位置表示界面。

优选的，对于多介质，需要满足不存在3种或3种以上有共同交接点的情况。

优选的，初始化流场qi,j,k的具体方法为：

根据不同介质的初始条件和介质m的标识初始化流场

优选的，区域表示介质m的虚拟流体之外使得计算区域规则的额外区域。

优选的，水平集方程为t表示时间，v表示速度矢量，表示梯度，φ表示符号距离函数。

优选的，更新介质m在n时刻的界面条件的方法为：

沿界面法向构造并求解多介质黎曼问题并求解：

un表示法向守恒变量，fn表示法向守恒通量，xn为界面法向位置，xi,n表示界面，求解获得

优选的，步骤(10)中n+1时刻整个计算域的流场为：

优选的，步骤(10)计算n+1时刻整个计算域的流场需要使用存储流场计算结果的数组q，存储n时刻全流域流场变量的临时数组qt1，存储n+1时刻获得的所有介质累加的流场变量的临时数组qt2。

优选的，整个计算域的流场随时间推进的存储变化过程如下：

10.1利用qt1临时存储时刻n的全流域流场变量令m＝1；

10.2根据qt1内存储的全流域流场变量定义介质m的界面条件，根据介质m的界面条件和流场更新n+1时刻介质m的流场变量并存储到q；

10.3将q内流场变量与qt2内流场变量累加并存储到qt2；

10.4判断m是否等于nmed，如果等于则进入步骤10.5，否则令m＝m+1，返回步骤10.2。

10.5将数组qt2内存储的累加结果赋值到数组q。

本发明与现有技术相比，其优点和创新性主要体现在以下几个方面：

(1)本发明的仿真流程设计方法能够实现多介质流体界面运动的仿真，适用于各类多介质流体界面演变问题的求解和模拟。

(2)本发明的仿真流程通用性强：多介质流体界面运动问题的解算设计流程可以直接应用于一般的cfd软件，使之在不影响该软件已有功能和计算效率的基础上，具备模拟多介质流体界面问题的能力；

(3)本发明代码修改量小：不改变原有程序架构，其非定常计算模块均可利用，也不需要对界面附近网格和数值格式进行特殊处理，原则上只需通过添加几个独立的计算模块和修改少量代码即可实现；

(4)本发明内存需求低：将界面看成一类特殊的边界，边界外虚拟流场与另一种介质真实流场的重合区域变量更新时内存需求低，适用于多种介质的耦合计算，不随介质数量的增加而增大。

附图说明

图1为多介质流体界面运动仿真流程图；

图2为不同介质计算区域示意图；

图3为模拟单介质流体非定常问题与多介质流体界面问题的模块对应关系；

图4为界面附近流场变量存储方式。

具体实施方式

本发明主要基于一般的cfd非定常计算系统，在能够进行单介质流场非定常计算的基础上，实现对多介质流体界面问题的模拟，见图1，cfd非定常计算流程为：

前处理过程(i)：

i-1.起始阶段，读取流场中涉及的介质总数nmed，依次输入每种介质的物性参数em，以及原始变量的初始条件控制参数cm等计算参数(m＝1,...,nmed)，获得不同介质特性的计算需求；原始变量例如可以为压力、密度、速度矢量等。

i-2.输入流场计算网格信息；

i-3.输入边界条件类型；

i-4输入初始界面信息包括界面形状、位置等特征；

i-5计算网格单元几何信息；

i-6.根据初始界面信息初始化符号距离函数i,j,k表示网格节点。即计算出每个网格节点到界面的符号距离，每个网格对应一种流体，符号为正值的节点代表一种流体，符号为负值的节点代表另一种流体，符号为零的位置表示界面；对于多种介质，需要满足不存在3种或3种以上有共同交接点的情况，则多种介质分别采用符号进行区域性的区分。

i-7.根据符号距离函数获得每种介质初始时刻占据的流场区域并确定网格节点与每种介质初场对应关系的标识

i-8.根据不同介质的初始条件和介质标识初始化流场

解算及结果输出过程(ii)：

ii-1.计算时间步长δt，δt根据需要选择，例如可以采用δx、δy、δz表示三个方向网格间距，c为声速，u、v、w为速度矢量的三个方向分量。

ii-2.根据输入边界条件类型对边界条件赋值，获得计算域边界条件；

ii-3.确定每种介质的计算区域表示介质m的界面外虚拟流体区域(一般跨过界面至少两个网格，与采用的数值格式有关)，表示介质m的虚拟流体之外使得计算区域规则的额外区域(保留求解空间，对计算结果无影响，便于每种介质流场计算)，见图2；

ii-4.求解水平集方程

获得新时刻符号距离函数将界面信息推进到新时刻，t表示时间，其中v表示速度矢量，表示梯度；

ii-5.根据新时刻的符号距离函数获得每种介质占据的流场区域并确定网格节点与每种介质对应关系的标识

ii-6.更新界面条件。即对虚拟流体区域中的每一个节点，沿界面法向构造并求解多介质黎曼问题

获得界面压力、速度和界面两侧密度,更新介质m在该节点的界面条件这里，u和f是守恒变量和通量，由qi,j,k采用现有公式计算获得；下标“n”表示界面法向，“-”和“+”分别表示界面附近位于φ<0和φ>0区域的流场状态；

ii-7.在更新界面条件后，根据每种介质的物性参数em和控制参数cm，在计算区域中依次对每种介质的守恒通量进行空间离散、对守恒变量进行时间离散，计算下一个时刻的守恒变量，根据下一个时刻整个计算域的守恒变量计算流场将每种介质的流场推进到新时刻；

ii-8.根据新时刻不同介质的流场和标识将各种介质的计算结果整合到整个计算域，获得新时刻流场

ii-9.判断是否到达结果输出时间，如果到达则进入ii-10，否则返回ii-1；

ii-10.输出新时刻流场计算结果通过可以显式出多介质流场特性；

ii-11.输出新时刻符号距离函数计算结果通过可以显示出多介质流体界面；

ii-12.判断时间推进是否结束，如果结束则完成多介质流体界面运动仿真；否则返回ii-1。

本发明将原始的非定常计算程序划分成2类模块：a类模块和b类模块。a类模块是通用模块，如边界条件输入模块、时间步长计算模块、结果输出模块等，在求解不同问题时共享使用，无需改动；b类模块需要对源代码进行少量修改，在保留原始功能的基础上满足界面问题解算需求。见图3。还需要添加一些c类模块，只有在计算多介质流体界面运动问题时，这些模块会被调用。对于已经具有a类模块和b类模块的非定常计算程序，针对b类和c类模块，设计模块修改和添加方式：

(1)b类模块：修改物性参数和计算参数输入模块(i-1)、流场初始化模块(i-8)、流场解算模块(ii-7)，在原始代码中增加次数为nmed的循环控制语句。当nmed为1时，表示只有一种介质，此时程序恢复成原始功能和流程。修改流场解算模块(ii-7)，依次传递每种介质的计算区域参数和状态方程参数；

(2)c类模块：添加初始界面信息输入模块(i-4)、符号距离函数初始化模块(i-6)、网格节点与介质流场对应关系标识模块(i-7、ii-5)、介质计算区域设定模块(ii-3)、符号距离函数解算模块(ii-4)、界面条件解算模块(ii-6)、多介质解算结果整合模块(ii-8)、符号距离函数输出模块(ii-11)。

本发明对界面附近变量存储的内存需求极小，而且不随介质数量的增加而增大。下面以一维问题为例具体说明，见图4。假设时刻n介质m和介质m+1的界面位于节点i-1和节点i之间。新时刻n+1界面变化到节点i和节点i+1之间。q是存储全流域流场变量的数组，在仿真过程中，也用于临时存储每个介质的新时刻流场变量仅需要两个临时数组qt1和qt2即可以完成多介质流场变量更新的操作，使q存储新时刻全流域流场变量

(1)利用qt1临时存储时刻n的全流域流场变量q：qt1←q；

(2)根据qt1定义介质m的界面条件，并更新介质m流场变量q到新时刻n+1，此时q内存储为介质m的流场变量

(3)将时刻n+1介质m流场变量q存储到qt2：qt2←q；

(4)令m＝m+1，重复步骤(2)-(3)直到所有介质流场变量都整合到qt2；

(5)更新时刻n+1的全流域流场变量q←qt2。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。