一种基于OpenFOAM的多层流体晃荡的模拟方法

一种基于openfoam的多层流体晃荡的模拟方法
技术领域
1.本发明属于多层流体晃荡技术领域，具体涉及一种基于openfoam的多层流体晃荡的模拟方法。


背景技术：

2.液舱晃荡问题属于共振问题的一种，几乎任何储存具有自由表面液体的移动结构都必须考虑晃荡问题，它是指在有限区域内两种或两种以上互不相溶的流体，由于外部激励接近液体共振频率而引起的有限区域内液体剧烈运动。
3.cfd方法将流体力学控制方程以离散化的代数形式表示，通过计算机求解离散化后的代数方程组，获得时空间上点的数值解。通过cfd方法，获取液舱的运动与受力情况，能够作为带有液舱的海洋结构物的设计参考。
4.近年来，越来越多的浮式生产储油卸油平台(fpso)被投入到海洋油气资源的开采当中。在海洋油气资源储卸油方式中，干式储油技术很大程度上限制了水下储油系统的有效储备量，因此油水置换湿式储油技术被广泛应用于油气开发平台上。但若将湿式储油技术应用于fpso平台，其技术难点之一就在于平台内部储存的多层流体的晃荡现象，这对整个储油系统的结构设计和工艺设计都是一个巨大的挑战，因此有必要对浮式储油平台中存在的气-油-水多层流体晃荡现象进行研究分析。
5.然而，目前的液舱晃荡模拟研究主要集中在单层液体，对于单层液体的晃荡模拟研究技术已较为成熟，但对于双层液体甚至三层液体的晃荡研究起步较晚，迄今，一些研究者已基于势流理论对多层流体晃荡的运动机理进行了研究，相关的晃荡模型试验也有序进行，但针对多层流体晃荡的模拟研究较少。


技术实现要素：

6.本发明的主要目的在于克服现有技术的缺点与不足，提出一种基于openfoam的多层流体晃荡的模拟方法。
7.为了达到上述目的，本发明采用以下技术方案：
8.一种基于openfoam的多层流体晃荡的模拟方法，包括以下步骤：
9.s1、建立三维液舱模型，对液舱进行网格划分，并设置液舱各面的几何边界条件；
10.s2、对液舱中多相流的流体域进行设置，定义各相流体的初始位置；
11.s3、构建用于模拟液舱晃荡的动网格模型，设置液舱六自由度运动属性；
12.s4、对液舱晃荡的各项计算控制参数与模型属性参数进行设定；
13.s5、计算程序运行求解并导出用于后处理的计算数据。
14.进一步的，步骤s1中，建立三维液舱模型具体为：
15.在openfoam中建立液舱模型，液舱由4个侧壁面、底面及顶面，共6个面围成；
16.设置构成液舱的8个节点坐标，定义液舱所处空间位置；
17.设置点位置矢量缩放因子，确定点位置数据单位，进一步确定液舱具体尺寸；
18.对液舱进行网格划分具体为，设置3个坐标轴垂直方向的网格数量及网格膨胀率，完成液舱内背景网格的设置。
19.进一步的，设置液舱各面的几何边界条件具体为：
20.除顶面外的其余5个面均设置为wall边界条件，对于顶面，若液舱顶部封闭，设置为wall边界条件，若液舱顶部开放，设置为patch边界条件。
21.进一步的，步骤s2具体为：
22.采用多层流体体积法实现多相流各流体域设定及多层液面同时捕捉；
23.设置未给定区域为空气相流体域，将空气相的相体积分数设置为1，其余相流体的相体积分数设置为0；
24.设置除空气相外的其余各相流体的流体域，确定多相流晃荡中各相流体所在位置；
25.将每一相流体所处流体域中，自身的相体积分数设置为1，其余相流体的相体积分数设置为0，完成各相流体的流体域的给定。
26.进一步的，步骤s3具体为：
27.基于openfoam单方向线性简谐运动动网格求解器程序，开发实现六个自由度方向同时运动的六自由度简谐运动动网格求解器程序；
28.通过动网格功能，施加外部简谐激励，实现液舱的刚体简谐运动；
29.设置液舱结构重心，并设置液舱六自由度方向运动的运动幅值及运动圆频率，实现动网格的运动方式构建。
30.进一步的，步骤s4中，需要设定的各项计算控制参数与模型属性参数具体包括：
31.重力大小及方向、流体属性、湍流模型、控制参数、并行参数、矩阵求解器以及离散格式。
32.进一步的，步骤s4具体为：
33.设置重力的大小及重力方向；
34.设置各相流体的具体流体属性，包括运动粘度、流体密度、表面张力系数以及是否为牛顿流体；
35.湍流模型设为层流假设、雷诺平均或大涡模拟；
36.根据选用的湍流模型不同，需设置不同的物理边界条件，层流假设需要设置初始速度场、初始压力场及各相流体相体积分数边界条件，雷诺平均需要在层流假设的物理边界条件基础上追加设置湍流脉动动能、湍流耗散率及湍流粘度边界条件，大涡模拟需要在层流假设的物理边界条件基础上追加设置湍流脉动动能、湍流粘度及模型中间变量边界条件；
37.将选用的cfd求解器设置为openfoam的multiphaseinterfoam求解器，并设置用于该求解器的控制参数，包括计算开始的物理时间、计算结束的物理时间、计算的时间步长、计算结果输出的时间间隔以及计算数据输出的格式；
38.设置后处理计算工具，采用probes探针工具记录壁面某点处压力场数值，采用基于interfaceheight界面高度工具开发的用于动网格界面高度记录的工具movinginterfaceheight记录各相流体的液面高度数据；
39.设置并行计算使用处理器数，并设定区域分解方法，完成并行运算的区域划分；
40.设置各物理量采用的离散矩阵求解器、ns方程求解算法及ns方程各项的离散方式。
41.进一步的，步骤s5具体为：
42.程序运行求解，输出用于后处理的数据结果；
43.后处理结果分为两个部分，一部分为postprocess后处理计算工具结果，包括液面高程及壁面压力；另一部分为模拟的场结果；
44.通过paraview实现后处理数据的可视化。
45.进一步的，开发实现六个自由度方向同时运动的六自由度简谐运动动网格求解器程序具体为：
46.在oscillatinglinearmotion动网格程序基础上开发harmonicmotion动网格程序，包括：
47.建立头文件harmonicmotion.h，在头文件中定义动网格六自由度运动所需的各项参数，包括结构重心矢量1个、六自由度运动幅值标量6个、六自由度运动圆频率标量6个，共13个参数；
48.建立定义文件harmonicmotion.c，在成员函数transformation()中构建网格运动函数，具体为：
49.首先提取当前计算时刻标量，然后通过当前时刻以及头文件中定义的运动幅值与圆频率，进行沿坐标轴方向简谐运动的计算，定义矢量记录平移运动，进行绕坐标轴方向简谐运动的计算，定义四元数记录旋转运动，最后定义七元数，结合头文件中定义的结构重心，整合网格的六自由度简谐运动，完成程序搭建。
50.本发明与现有技术相比，具有如下优点和有益效果：
51.1、本发明方法建模和计算工作都采用开源软件openfoam完成，可根据模拟需要对源程序代码进行直接修改和补充，避免了采用商业软件的购置成本高、无法进行二次开发等缺点。
52.2、本发明方法在用于捕捉自由液面的vof方法的基础上，采用多层vof方法，实现了多层流体的模拟及多液面的同时捕捉，从而实现多层流体晃荡的模拟。
53.3、本发明方法克服了oscillatinglinearmotion动网格程序仅能实现单自由度的简谐运动的束缚，开发了能够同时实现六自由度简谐运动的harmonicmotion动网格程序，能够实现多方向液舱晃荡的模拟。
附图说明
54.图1是本发明方法的流程图；
55.图2是多层流体晃荡模型示意图；
56.图3是液舱网格划分情况示意图；
57.图4是液舱中水相流体域示意图；
58.图5是液舱中油相流体域示意图；
59.图6是多层流体晃荡的波高时历曲线图；
60.图7是多层流体晃荡的壁面压力时历曲线图；
61.图8是多层流体晃荡4s时刻波面图；
62.图9是多层流体晃荡4.5s时刻波面图；
63.附图标号说明：1-液舱；2-空气相；3-油相；4-水相。
具体实施方式
64.下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
65.实施例
66.本实施例基于les湍流模型模拟液舱中含自由表面的双层液体，在外部水平简谐激励下的三维晃荡运动，激励频率接近上层液体的共振频率。液舱尺寸为0.5m
×
0.1m
×
0.5m，下层水和上层油高度均为0.075m，外部激励圆频率为7.3263rad/s(上层液体共振频率的1.1倍)，激励幅值为0.001m，运动方程为x(t)＝0.001sin7.3263t。
67.如图1所示，本发明，一种基于openfoam的多层流体晃荡的模拟方法，包括以下步骤：
68.s1、建立三维液舱模型，对液舱1进行网格划分，并设置液舱各面的几何边界条件；具体包括：
69.s11、在openfoam中建立液舱模型，液舱由4个侧壁面、底面及顶面共6个面围成。
70.s12、设置构成液舱的8个节点坐标，定义液舱所处空间位置。在vertices中将0～7号节点分别设置为(-0.25
ꢀ‑
0.05 0)，(-0.25 0.05 0)，(-0.25 0.05 0.5)，(-0.25
ꢀ‑
0.05 0.5)，(0.25
ꢀ‑
0.05 0)，(0.25 0.05 0)，(0.25 0.05 0.5)，(0.25
ꢀ‑
0.05 0.5)；
71.s13、设置点位置矢量缩放因子，确定点位置数据单位，进一步确定液舱具体尺寸；在convertometers中将点位置矢量缩放因子设置为1，即设定单位为m；
72.s14、设置3个坐标轴垂直方向的网格数量及网格膨胀率，完成液舱内背景网格的设置；在本实施例中，具体为：
73.在blocks中将网格数量与网格膨胀率设置为hex(0 1 2 3 4 5 6 7)(20 125 200)simplegrading(1((0.5 0.8 1)(0.5 0.2 1))1)；即将由0～7号节点构成的块划分网格，宽度方向网格数20，高度方向网格数125，长度方向网格数200，宽度与长度网格膨胀率为1，即均匀网格，高度方向网格非均匀，其中靠下50％的部分占有80％的网格，靠上50％的部分占有20％的网格。
74.s15、设置三维液舱6个面的几何边界条件，除顶面外的其余5个面均设置为wall边界条件，对于顶面，若液舱顶部封闭，设置为wall边界条件，若液舱顶部开放，设置为patch边界条件；在本实施例中，具体为：
75.在boundary中将几何边界条件设置为top(2 3 7 6)wall，bottom(0 1 5 4)wall，left(0 1 2 3)wall，right(4 5 6 7)wall，front(0 4 7 3)wall，back(1 5 6 2)wall。将液舱的6个面均设置为可识别壁面边界条件wall。
76.打开paraview可视化后处理软件，可在paraview中查看网格情况，网格划分情况如图3所示，网格数为500000。
77.s2、对液舱中多相流的流体域进行设置，定义各相流体的初始位置；具体为：
78.s21、采用多层流体体积法实现多相流各流体域设定及多层液面同时捕捉。流体体积法通过互不相容的流体组分共用一套动量方程，引入相体积分数表示某一相的体积占所
在网格体积的比值从而实现两相流的模拟，而多层流体体积法在此基础上通过引入多个相体积分数，实现多相流的模拟。具体使用时，选用openfoam的multiphaseinterfoam求解器来实现。
79.s22、设置未给定区域为空气相流体域，将空气相2的相体积分数设置为1，其余相流体的相体积分数设置为0。
80.在本实施例中，具体为：
81.在defaultfieldvalues中将volscalarfieldvalue设置为alpha.air 1、alpha.water0以及alpha.oil 0。即默认区域为空气相，该相中空气的相体积分数为1，水的相体积分数为0，油的相体积分数为0。
82.s23、设置除空气相外的其余各相流体的流体域，确定多相流晃荡中各相流体所在位置。
83.在本实施例中，具体为：
84.在regions中的boxtocell中将box设置为water(-1
ꢀ‑
1 0)(1 1 0.075)，oil(-1
ꢀ‑
1 0.075)(1 1 0.15)。即水相流体域为(-1
ꢀ‑
1 0)和(1 1 0.075)两个点构成的体对角线所确定的长方体区域，油相流体域为(-1
ꢀ‑
1 0.075)和(1 1 0.15)两个点构成的体对角线所确定的长方体区域。
85.s24、将每一相流体所处流体域中，自身的相体积分数设置为1，其余相流体的相体积分数设置为0，完成各相流体流体域的给定。将水相4的volscalarfieldvalue设置为alpha.air0、alpha.water1、alpha.oil 0，将油相3的volscalarfieldvalue设置为alpha.air0、alpha.water0、alpha.oil1。
86.至此，完成多相流设置工作，打开paraview可视化后处理软件，在paraview中查看多相流设置情况，多相流设置情况如图4和图5所示。
87.如图2所示，是多层流体晃荡模型示意图。
88.s3、构建用于模拟液舱晃荡的动网格模型，设置液舱六自由度运动属性；具体为：
89.s31、基于openfoam单方向线性简谐运动动网格求解器程序，开发六自由度简谐运动动网格求解器程序。在oscillatinglinearmotion动网格程序基础上开发harmonicmotion动网格程序，在本实施例中，具体操作为：
90.建立头文件harmonicmotion.h，在头文件中定义动网格六自由度运动所需的各项参数，包括结构重心矢量1个、六自由度运动幅值标量6个、六自由度运动圆频率标量6个，共13个参数。
91.建立定义文件harmonicmotion.c，在成员函数transformation()中构建网格运动函数，首先提取当前计算时刻标量，然后通过当前时刻以及头文件中定义的运动幅值与圆频率，进行沿坐标轴方向简谐运动的计算，定义矢量记录平移运动，同样的，进行绕坐标轴方向简谐运动的计算，定义四元数记录旋转运动，最后定义七元数，结合头文件中定义的结构重心，整合网格的六自由度简谐运动，完成程序搭建。
92.完成程序搭建后，通过wmake进行编译，最终程序能够网格在坐标轴六个自由度方向上同时进行简谐运动，相较软件自带的单方向简谐运动更具实用性。
93.至此，完成harmonicmotion动网格简谐运动程序的构建。
94.s32、通过动网格功能，施加外部简谐激励，实现液舱的刚体简谐运动。
95.在本实施例中，将dynamicfvmesh设置为dynamicmotionsolverfvmesh，将motionslover设置为solidbody，将solidbodymotionfunction设置为harmonicmotion。
96.s33、设置液舱结构重心，并设置液舱六自由度方向运动的运动幅值及运动圆频率，实现动网格的运动方式构建。
97.在本实施例中，即液舱做纵荡简谐运动，运动幅值为0.001m，运动圆频率为7.3263rad/s。
98.s4、对液舱晃荡的各项计算控制参数与模型属性参数进行设定，包括重力方向、流体属性、湍流模型、控制参数、并行参数、矩阵求解器、离散格式；具体为：
99.s41、设置重力的大小及重力方向，并设置各相流体的具体流体属性，包括运动粘度、流体密度、表面张力系数以及是否为牛顿流体，完成物理参数的设置。
100.在本实施例中，将重力g的大小设置为9.81m/s2，方向为z轴负方向；将水相设置为牛顿流体，运动粘度为1
×
10-6
m2/s，流体密度为1000kg/m3；油相设置为牛顿流体，运动粘度为5
×
10-6
m2/s，流体密度为846kg/m3；空气相设置为牛顿流体，运动粘度为1.48
×
10-5
m2/s，流体密度为1kg/m3。将各相流体间的表面张力系数sigma均设置为默认值0.07。
101.s42、设置选用的湍流模型，可以为层流假设、雷诺平均或大涡模拟。
102.在本实施例中，采用大涡模拟方法进行湍流模拟，将model设置为keqn，使用一方程涡粘模型作为les湍流模型，将turbulence设置为on，启用湍流模型，将printcoeffs设置为on，输出湍流模型系数，将delta设置为smooth，代表delta模型选用光顺模型。
103.s43、根据选用的湍流模型不同，需要设置不同的物理边界条件，大涡模拟需要初始速度场、初始压力场、各相流体相体积分数、湍流脉动动能、湍流粘度及模型中间变量边界条件。
104.在本实施例中，将初始速度场u的边界条件设置为移动壁面速度边界条件movingwallvelocity，初始速度边界值为0m/s；将初始压力场p_rgh的边界条件设置为固定压力通量边界条件fixedfluxpressure，初始压力边界值为0kg/(m
·
s2)；将水相、油相、空气相的相体积分数alpha的边界条件设置为零法相梯度边界条件zerogradient，初始相体积分数为0；将湍流脉动动能k的边界条件设置为给定数值边界条件fixedvalue，初始湍流脉动动能边界值为1
×
10-11
m2/s2；将湍流粘度nut的边界条件设置为零法相梯度边界条件zerogradient，初始湍流粘度边界值为0m2/s；将模型中间变量nutilda的边界条件设置为给定数值边界条件fixedvalue，初始中间变量边界值为0m2/s。
105.s44、将选用的cfd求解器设置为openfoam的multiphaseinterfoam求解器，并设置用于该求解器的控制参数，包括计算开始的物理时间，计算结束的物理时间，计算的时间步长，计算结果输出的时间间隔，计算数据输出的格式。
106.在本实施例中，选用multiphaseinterfoam多相流求解器。于物理时间0s时刻开始计算，于物理时间100s时刻结束计算，时间步长为0.005s。每0.25s输出一次场结果，以ascii编码形式输出结果。
107.s45、设置案例采用的后处理计算工具，采用probes探针工具记录壁面某点处压力场数值，采用基于interfaceheight界面高度工具开发的用于动网格界面高度记录的工具movinginterfaceheight记录各相流体的液面高度数据。继续在controldict文件中进行设置，在functions中加入probe，将probe中的type设置为probes，使用探针工具。将
writecontrol设置为timestep，表示每个时间步都输出一个结果。将probelocations设置为((0.25 0 0.075)(0.25 0 0.15))，表示探针坐标为位于壁面上的点。将fields设置为(p)，表示探针用于记录压力场数据。
108.基于interfaceheight工具，构建能够实现测点与动网格同步运动的波高记录工具movinginterfaceheight，在本实施例中，具体操作为：
109.建立头文件movinginterfaceheight.h，其中无需新增参数定义，只需将interfaceheight关键字均改为movinginterfaceheight即可。
110.建立定义文件movinginterfaceheight.c，在成员函数writepositions()中加入测点运动函数，首先提取当前计算时刻，之后根据所需实现的简谐运动规律，对测点坐标三维数组的值进行修改，使其与网格同步运动
111.完成程序搭建后，通过wmake进行编译，最终程序能够实现设置测点与动网格以同样的规律进行运动，保证测点相对于随体坐标是不动的，保证测点记录数值的准确性。
112.至此，完成movinginterfaceheight波高监测程序的构建。
113.完成程序构建后，回到controldict中，在functions中加入waterinterfaceheight与airinterfaceheight，分别用于记录下层液体波高与上层液体波高。将type均设置为movinginterfaceheight，表示使用自主开发的移动波高记录工具。将alpha分别设置为alpha.water与alpha.air，表示通过相体积分数为0.5时确认波面位置，从而记录某相的高度，由于oil相有两处alpha为0.5的位置，因此使用air相间接记录oil相的自由表面位置。将locations设置为(0.25 0 0.2)，表示波高测点位置，如此，即可使得波高记录测点相对随体坐标而言保持在(0.25 0 0.2)的位置，从而记录液舱壁面处的波高变化情况。值得注意的是，movinginterfaceheight记录的是液体表面处距离壁面的距离，waterinterfaceheight为下层液体表面到下壁面的距离，airinterfaceheight为上层液体表面到上壁面的距离，需要另外经过换算才正式得到波高。
114.s46、设置并行计算使用处理器数，并设定区域分解方法，完成并行运算的区域划分。
115.本实施例中，使用16个处理器进行并行运算。采用scotch方法完成区域分解。
116.s47、设置各物理量采用的离散矩阵求解器、ns方程求解算法及ns方程各项的离散方式。
117.在本实施例中，fvsolution文件中，solvers为各物理量分别采用的矩阵求解器，p_rgh采用gamg求解器，p_rghfinal采用pcg求解器，u|k采用smoothsolver求解器。tolerance与reltol分别为残差与相对残差。preconditioner为预条件处理器。pimple下为采用的pimple算法参数。relaxationfactors下为松弛参数设置，松弛因子设定为1，即不松弛。
118.设置ns方程各项的离散格式。ddtschemes为时间导数项，采用euler格式，gradschemes为压力梯度项，采用gausslinear格式，divschemes为对流项，其中的各项均有各自的离散格式，laplacianschemes为扩散项，采用gausslinearcorrected格式；而interpolationschemes为网格体心插值到面心所用离散格式，采用linear格式；sngradschemes为面的法向相关的离散格式，采用corrected格式。
119.s5、计算程序运行求解并导出用于后处理的计算数据；具体包括：
120.s51、程序运行求解，输出用于后处理的数据结果。
121.s52、后处理结果分为两个部分，一部分为postprocess后处理计算工具结果，包括液面高程及壁面压力。将postprocessing文件夹中的波高与压力结果绘制成曲线，如图6和图7所示；另一部分为模拟的场结果。
122.通过paraview实现后处理数据的可视化，如图8和图9所示。
123.还需要说明的是，在本说明书中，诸如术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
124.对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其他实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。