一种三维多相可压缩流固耦合的快速计算方法与流程

1.本发明属于流体计算技术领域，具体涉及一种三维多相可压缩流固耦合的快速计算方法。


背景技术：

2.由于水下强冲击是个多相、多介质、多物理场的高度复杂跨学科物理化学问题，模拟该类问题采用的技术包括高阶精度激波捕捉算法、可压缩多相计算流体力学算法、大变形大位移结构有限元算法、强冲击多介质流固耦合算法。目前在解决水下强冲击方面国际上并没有特别高效与稳定的算法。虽然研究人员在局部领域都取得了一定的突破，但是将这些算法完全耦合在一起形成一套完整高效的水下冲击模拟算法还存在很多问题，包括多相流界面解耦、流固界面解耦、流固界面空化塌陷，以及激波与气泡动力学阶段整体模拟的计算量巨大的问题。
3.目前应用于水下冲击的主要方法包括任意欧拉-拉格朗日耦合算法，浸入边界方法以及松流固耦合算法。
4.任意欧拉-拉格朗日耦合算法通用性好，是解决水下冲击可压缩流固耦合目前应用性比较广的方法，但是该方法存在流固界面处理复杂，在结构变形的情况下需要重构网格导致计算效率和稳定性下降的问题。
5.浸入边界方法具有流固界面处理简单以及不需要重构网格的优势。目前浸入边界方法主要应用于不可压缩流体与结构耦合的问题求解，在水下冲击的多相可压缩流固耦合领域的应用还存在多相流动与可压缩流动求解的技术难点。


技术实现要素：

6.本发明的目的在于提供一种三维多相可压缩流固耦合的快速计算方法，克服了现有技术的不足，通过开发基于虚拟介质方法的三维多相可压缩流体算法与流固耦合算法大幅度提高了水下冲击模拟计算的精度与稳定性；开发水下冲击激波阶段与气泡动力学阶段相互转换的算法，基于流场的压力对两种算法的转换进行判断，大幅度提高了水下冲击问题整体求解的效率。
7.为解决上述问题，本发明所采取的技术方案如下：
8.一种三维多相可压缩流固耦合的快速计算方法，包括以下步骤：
9.(1)建立受冲击材料的三维模型并划分笛卡尔非均匀网格；
10.(2)建立无粘可压缩流体在笛卡尔坐标下的质量、动量与能量方程以及状态方程；
11.(3)进行多相可压缩流固耦合计算；
12.(4)通过流固耦合求解器提供初始压力给结构场，同时获得结构场提供的初始节点速度；
13.(5)通过可压缩-不可压缩转换求解算法对水下冲击全周期和全流程进行模拟，输出计算结果。
14.进一步，所述步骤(2)中无粘可压缩流体在笛卡尔坐标下的质量、动量与能量方程，可表示为：
[0015][0016]
其中，
[0017]
ρ是流体密度,p是压力,u,v,w是x,y,z方向的速度,g重力加速度；
[0018]
总能量e
t
＝e 0.5(u2 v2 w2)，其中，e是内能。
[0019]
进一步，所述步骤(2)中状态方程可表示为
[0020]
p＝p(ρ,e),ρ＝ρ(p,e),ore＝e(ρ,p).。
[0021]
进一步，所述步骤(3)中多相可压缩流固耦合计算，具体包括
[0022]
通过距离函数确定界面的位置；
[0023]
定义虚拟节点的虚拟流体状态以及在界面进行近似黎曼的求解。
[0024]
进一步，所述距离函数可以通过求解如下瞬态运输方程获得：
[0025]
φ
t
uφ
x
vφy wφz＝0
[0026]
ψ
t
uψ
x
vψy wψz＝0
[0027]
其中，u,v,w是流体在x,y,z方向的速度；φ
t
，φ
x
，φy，φz是第一个两相流界面的距离函数随时间与空间的变化量，ψ
t
，ψ
x
，ψy，ψz是第二个两相流界面的距离函数随时间与空间的变化量。
[0028]
进一步，所述近似黎曼的求解可表示为：
[0029][0030][0031][0032][0033]
其中，下标“i”,“il”和“ir”指的界面，以及界面左右两侧；ρ
il
(ρ
ir
)与c
il
(c
ir
)是界面左右两侧的流体密度与声速；ui，pi为界面的速度与压力；u
il
，u
ir
可以通过界面的两条可回归到相应介质的非线性特征线进行插值获得，或者设置为u
il
＝u
i-1
，u
ir
＝u
i 2
；ρ
l
(pi)，ρr(pi)是在界面压力pi下左右激波流动的流体密度。
[0034]
进一步，所述步骤(4)中流固耦合求解器为耦合欧拉函数和拉格朗日函数的求解器。
[0035]
进一步，所述步骤(5)中可压缩-不可压缩转换求解算法的流程包括：在转换时，可压缩算法在所有体积网格里都有流动和流体的物性参数；采用界面跟踪算法计算界面的位置，包括气泡与自由界面的位置；边界元方法根据可压缩流体计算的结果产生理想气泡以及自由界面的位置；产生势流边界元网格；将转换时可压缩流体的流场进行插值获得边界元网格法向的速度场；应用上述流程在流固界面上；格林方程应用于求解在bem网格上的势；进行bem的计算，并输出结果。
[0036]
本发明与现有技术相比较，具有以下有益效果：
[0037]
1、本发明通过耦合基于虚拟介质方法与高阶激波捕捉方法对物质界面进行高精度的解耦，对物质界面波的结构进行精确的求解；确保三相多相可压缩流固耦合算法的精度和稳定性。
[0038]
2、本发明通过水下冲击激波阶段与气泡动力学阶段相互转换的算法，基于流场的压力对两种算法的转换进行判断，大幅度提高了水下冲击问题整体求解的效率。
附图说明
[0039]
图1为多相流界面示意图。
[0040]
图2为虚拟流体节点示意图。
[0041]
图3为边界逼近法搜索流体网格到节点距离、流体网格到线段距离、流体网格到平面的最短距离的示意图。
[0042]
图4为流固耦合界面的示意图。
[0043]
图5为流体网格速度的示意图。
[0044]
图6为受冲击复合材料复合材料螺旋桨的几何拓扑结构。
[0045]
图7为受冲击复合材料复合材料螺旋桨的计算网格。
[0046]
图8为激波阶段不同时间下流场的压力云图以及螺旋桨的变形。
[0047]
图9为气泡阶段不同时间下气泡的演变过程以及螺旋桨的变形。
具体实施方式
[0048]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0049]
本发明所述一种三维多相可压缩流固耦合的快速计算方法，包括以下步骤：
[0050]
(1)建立受冲击材料的三维模型并划分笛卡尔非均匀网格；
[0051]
(2)建立无粘可压缩流体在笛卡尔坐标下的质量、动量与能量方程以及状态方程；
[0052]
(3)进行多相可压缩流固耦合计算；
[0053]
(4)通过流固耦合求解器提供初始压力给结构场，同时获得结构场提供的初始节点速度；
[0054]
(5)通过可压缩-不可压缩转换求解算法对水下冲击全周期和全流程进行模拟，输出计算结果。
[0055]
1、无粘可压缩流体在笛卡尔坐标下的质量、动量与能量方程，可表示为：
[0056][0057]
其中，
[0058]
ρ是流体密度,p是压力,u,v,w是x,y,z方向的速度,g重力加速度；
[0059]
总能量e
t
＝e 0.5(u2 v2 w2)，其中，e是内能。
[0060]
2、状态方程可表示为：
[0061]
p＝p(ρ,e),ρ＝ρ(p,e),ore＝e(ρ,p).。
[0062]
通常情况下，状态方程可以用mie-gruneisen的统一格式进行描述，各种介质的状态方程可以如下表所示：
[0063][0064]
3、如图1所示，两个界面分割了三种介质，第一个levelset函数(φ＝0)定义了介质1(med1)和介质2(med2)之间的界面，同时第二个levelset函数(ψ＝0)定义了介质2(med2)和介质3(med3)之间的界面。
[0065]
每一个介质的定义如下所示：
[0066][0067]
这里med2通常被视为主介质，med1&med2组成了第一个两相流，而med3&med2组成
了另外一个两相流。由于med1和med2的计算相互独立，因此这些两种两相流解的组合将给出三相流动最终的解。对于每种介质，距离函数可以通过求解如下瞬态运输方程获得：
[0068]
φ
t
uφ
x
vφy wφz＝0
[0069]
ψ
t
uψ
x
vψy wψz＝0
[0070]
其中，u,v,w是流体在x,y,z方向的速度；φ
t
，φ
x
，φy，φz是第一个两相流界面的距离函数随时间与空间的变化量，ψ
t
，ψ
x
，ψy，ψz是第二个两相流界面的距离函数随时间与空间的变化量。
[0071]
如果两个界面重合，这计算med1与med2就是相关的，在这个场景下，计算med1/med3需要来自于med3/med1的信息，因此本专利方法通过以下方程来定义不同的介质：
[0072][0073]
4、如图2所示，假设材料界面在时间t＝tn位于网格节点i和i 1之间，我们要计算下一个时间t＝t
n 1
的流场。界面的新位置首先使用距离函数的levelset技术获得，为了计算med1的流场，根据所使用算法的精度，网格节点i 1、i 2或更多节点需要被定义为med1虚拟流体节点。虚拟流体节点的状态通过与med1相同的状态方程来定义，应用相同的方法可以定义med2的虚拟流体节点(i,i-1,i-2...)和虚拟流体节点状。
[0074]
虚拟节点的虚拟流体状态可以通过在两相流界面的两个非线性特征方程在界面处定义和求解一个多介质黎曼问题：
[0075]
along
[0076]
along
[0077]
这里下标“i”,“il”and“ir”指的界面，以及界面左右两侧；ρ
il
(ρ
ir
)与c
il
(c
ir
)是界面左右两侧的流体密度与声速；ui，pi为界面的速度与压力；u
il
，u
ir
可以通过界面的两条可回归到相应介质的非线性特征线进行插值获得，或者简单的设置为u
il
＝u
i-1
，u
ir
＝u
i 2
。由于当激波对界面进行冲击时，,u
il
(p
il
)与u
ir
(p
ir
)在穿过界面的时候是不连续的(间断)，因此必须在界面进行近似黎曼的求解，如下所示：
[0078][0079][0080][0081]
[0082]
ρ
l
(pi)，ρr(pi)是在界面压力pi下左右激波流动的流体密度。
[0083]
5、在流固耦合求解器中，第一步是确定距离函数，它实际上是从流体网格到结构界面的最短距离。流固耦合在每个时间步长的距离函数都是通过边界逼近方法计算。如图3所示，边界逼近法搜索流体网格到节点距离、流体网格到线段距离、流体网格到平面的最短距离。为了降低计算成本，距离搜索仅在结构表面的狭窄区域内进行。
[0084]
使用上面提供的距离函数，第二步是耦合欧拉函数和拉格朗日函数的求解器。流体计算是在固定的笛卡尔网格上进行的，而结构计算是使用与结构有限元计算完成的。如图4所示流固耦合界面最终是形成阶梯近似。
[0085]
对于流固耦合求解器的关键是定义结构内部的虚拟单元，以便使用虚拟节点坐标和速度进行流体计算。每个流体网格的法线方向可以基于距离函数进行计算。结构内部的虚拟节点可以通过以下方程进行投影：
[0086][0087]
其中，q＝(ρ,p,u,v,w)是外插量的数组。沿法线方向采用简单的迎风格式计算τ，当达到稳态时，在虚拟节点的非物理的外插速度场可以通过以下方程进行重构：
[0088][0089]
其中是通过实际流体网格速度通过外插的方法获得的速度，vs是固体界面速度，如图5所示，该速度是在交叉点上的速度。
[0090]
6、可压缩算法可以模拟包括冲击波和气泡动力学，然而这种方法的计算成本和cpu时间导致无法在合理时间内对水下冲击全周期和全流程进行模拟。本专利开发一种组合转换求解算法，使用可压缩流体算法来模拟爆炸的初始阶段，然后切换到边界元方法来模拟气泡动力学，在气泡破裂时再切换回到可压缩流体算法。该转换求解算法的关键是转换时构建流场初始条件。
[0091]
以可压缩流动转换到不可压缩流动来说明这个算法的流程：
[0092]
(1)在转换时，可压缩算法在所有体积网格里都有流动和流体的物性参数；
[0093]
(2)采用界面跟踪算法计算界面的位置，包括气泡与自由界面的位置；
[0094]
(3)边界元方法根据可压缩流体计算的结果产生理想气泡以及自由界面的位置；
[0095]
(4)产生势流边界元网格；
[0096]
(5)将转换时可压缩流体的流场进行插值获得边界元网格法向的速度场；
[0097]
(6)应用(1)-(5)流程在流固界面上
[0098]
(7)格林方程应用于求解在bem网格上的势；
[0099]
(8)进行bem的计算，并输出结果。
[0100]
实施案例：
[0101]
(1)复合材料水下冲击案例-案例物理示意模型
[0102]
受冲击复合材料复合材料螺旋桨的几何拓扑结构如图6所示，该螺旋桨是本案例主要研究在受冲击下响应的结构。
[0103]
(2)复合材料水下冲击案例-计算网格
[0104]
复合材料受水下冲击案例的计算网格如图7所示，生成的是笛卡尔非均匀网格。为了精准的计算激波阶段，在冲击源物理区域生成较密的结构网格。
[0105]
(3)复合材料水下冲击案例-激波阶段计算
[0106]
激波阶段，主要是应用多相可压缩流固耦合计算方法，图8展示了在激波阶段不同时间下流场的压力云图以及螺旋桨的变形。
[0107]
(4)复合材料水下冲击案例-气泡动力学计算
[0108]
气泡动力学阶段，采用了可压缩-不可压缩转换算法，主要是应用边界元流固耦合计算方法，图9展示了在气泡阶段不同时间下气泡的演变过程以及螺旋桨的变形。应用该转换算法将计算时长从30天下降到4小时，效率提高两个数量级以上。
[0109]
对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。