一种应用于异构并行两相数值水池系统的全解耦投影方法

本发明属于计算流体力学领域，具体涉及一种应用于异构并行两相数值水池系统的全解耦投影方法。

背景技术：

1、数值水池是以船舶与海洋结构物水动力学性能研究为特定应用领域的虚拟仿真试验系统。数值水池的基本原理是利用计算机模拟流体流动，通过求解流体运动控制方程，模拟船舶与海洋结构物的运动和受力。通过用计算机软件实现甚至超越物理水池的功能，从而设计出具有优良性能的船舶与海洋结构物。

2、计算流体力学作为数值水池技术的核心组成部分，是一门运用计算机科学、数值方法和流体动力学等多学科知识，通过在计算机上进行数值模拟方法解决流体流动问题的学科。旨在准确、高效地捕捉流体行为，在航空航天、海洋工程、化工和医疗等领域都有着广泛而重要的应用。现有计算流体力学方程主要基于描述流体运动的不可压navier-stokes方程，其包括连续性方程(1)和动量方程(2)：

3、

4、

5、只需要求解这两个方程组就可以得到速度场和压力场的完整解。目前针对不可压navier-stokes方程的数值求解算法有很多，常见的有simple算法，piso算法，涡量流函数法、人工压缩方法以及投影算法等；但他们大多都需要求解大型线性方程组，这严重影响了计算效率。为此学者开始研究针对不可压流体的更简便的求解算法，例如基于压力修正的投影算法，其求解不可压navier-stokes方程时，往往将不可压navier-stokes方程解耦为如下3个方程：

6、中间速度方程：

7、当前压力方程：

8、当前速度方程：

9、中间速度方程和当前速度方程都能显式处理，但当前压力方程不含压力瞬态项，只能隐式计算。当采用有限体积法进行离散求解时，隐式格式的每个时间步都需要重新构造系数矩阵并求解，这往往是多相流数值计算效率瓶颈的主要原因。为了提高多相流数值模拟的计算效率，学者们开展了许多研究工作，主要可以分为下述3类方法：第1类是常系数矩阵时间步进算法求解navier-stokes和cahn-hilliard方程，避免每个时间步的系数矩阵构造，并结合预处理技术，从而实现大密度比两相流的高效计算；也可以通过引入一组辅助标量，将系数矩阵转变为定常矩阵，实现高效且无条件能量稳定的梯度流求解。第2类是加速线性方程组本身的收敛过程，采用多重网格法在粗网格和细网格之间反复转换，有效地提高了迭代算法的收敛速度；第3类是并行算法的设计，比如区域分解算法，主要思想是将计算域分解成多个区域并行计算。除了上述的3类方法外，自适应网格，和人工神经网络等技术也常被用于改善不可压多相流的计算效率。但这些算法大多需要隐式或半隐式求解不可压多相流，线性方程组的存在既占用了大量内存，也严重影响了计算效率，这个难题在三维模型数值计算中将更加突出。

技术实现思路

1、为解决上述问题，实现基于cpu+gpu程序架构加速的不可压多相流的并行计算，有效提升数值水池的计算效率，本发明提供了一种应用于异构并行两相数值水池系统的全解耦投影方法，包括采用三维建模软件构建数值水池系统模型；基于数值水池系统模型，采用全解耦投影方法进行数值模拟计算，输出计算结果；其中，采用全解耦投影方法进行数值模拟计算的方法为：基于cpu/gpu异构并行架构，采用显式算法对不可压navier-stokes方程进行求解，得到异构并行两相数值水池系统模型的数值模拟结果。

2、进一步的，所述cpu/gpu异构并行架构包括一个cpu和多个gpu；其中cpu负责管理数据在gpu间的传输，gpu负责每一时间步的数值计算。

3、进一步的，基于cpu/gpu异构并行架构的求解流程包括以下步骤：

4、s1.在cpu上进行初始化：设置计算域和最大时间步istep_max；对数值水池系统模型进行网格划分，设置两相物性参数，初始化u矩阵、p矩阵和相界面；

5、s2.在gpu上判断当前时间步istep是否小于最大时间步istep_max，若是，则执行步骤s3，若不是，则结束计算并执行步骤s6；

6、s3.在gpu上根据两相界面运动场景建立各流体质点的cahn-hilliard方程，采用二阶runge-kutta算法求解cahn-hilliard方程，得到各流体质点的重构密度和重构黏度；

7、s4.在gpu上基于重构密度和重构黏度建立各流体质点的navier-stokes方程，采用全解耦显式算法求解navier-stokes方程，得到各流体质点的速度场和压力场；

8、s5.在gpu上令istep＝istep+1，并返回步骤s2；。

9、s6.gpu向cpu发送计算数据并释放gpu显存，cpu对计算数据转换为tecplot数据格式输出并进行可视化展示。

10、进一步的，步骤s3在gpu上采用多相流界面捕捉算法求解任一流体质点的cahn-hilliard方程，得到重构密度和重构黏度，包括：

11、s31.在gpu上根据两相界面运动场景建立流体质点的cahn-hilliard方程，表示为：

12、

13、

14、其中，ψ表示化学势，φ表示相变量，和表示散度算子和拉普拉斯算子，ε表示界面宽度，m表示迁移率，u表示速度矢量，

15、s32.设置化学势ψ和相变量φ的边界条件为：

16、

17、

18、s33.基于边界条件求解cahn-hilliard方程，得到相变量φ的数值大小；

19、s34.根据求解出的相变量φ进行线性插值计算，得到重构的密度和动力黏度，表示为：

20、

21、

22、其中，ρ表示重构密度，ρ1表示纯气相密度，ρ2表示纯液相密度，μ表示重构黏度，μ1表示纯气相黏度，μ2表示纯液相黏度。

23、进一步的，步骤s4在gpu上采用全解耦显式算法求解任一流体质点的navier-stokes方程，得到流体质点速度场和压力场，包括：

24、s41.在gpu上基于重构密度和重构黏度建立各流体质点的navier-stokes方程，表示为：

25、

26、

27、其中，m表示质量通量，u表示速度矢量，ρ表示重构密度，p表示压力，f表示重力和表面张力，μ表示重构黏度，和表示散度算子和梯度算子；

28、s42.将navier-stokes方程解耦为如下的3个方程：

29、中间速度方程：

30、当前压力方程：

31、当前速度方程：

32、其中，fext表示外力项；

33、s43.基于中间速度方程显式求解得到中间速度u*；

34、s44.将中间速度u*代入当前压力方程后进行修正，采用一阶欧拉格式显式推进修正后的当前压力方程，最后显式并行求解得到当前压力pn+1；

35、s45.将当前压力pn+1代入当前速度方程显式求解得到当前速度un+1。

36、进一步的，步骤s44将中间速度u*代入当前压力方程后进行修正，表示为：

37、

38、采用一阶欧拉格式显式推进修正后的当前压力方程，循环求解k次：

39、

40、

41、其中，k表示当前压力方程循环求解次数，n表示当前所处的时间步。

42、本发明的有益效果：

43、本发明提出了求解navier-stokes和cahn–hilliard方程的并行有限体积投影算法，通过引入人工的压力瞬态项，对泊松方程进行不完全迭代，实现了不可压navier-stokes方程的全解耦显式计算。此外，结合异构并行的程序架构，本算法完全避免了线性方程组的求解，并且可以对每个网格点分别进行计算，具备极佳的并行性能。