一种适用于模拟水弹性问题的耦合数值计算方法技术

本发明专利技术适用于计算力学领域，尤其是涉及一种适用于模拟水弹性问题的耦合数值计算方法，本发明专利技术采用SPH方法用来模拟流场，发挥粒子类方法在模拟涉及自由交界面及界面变形问题上的优势；采用FEM方法来模拟弹性体，发挥网格类方法处理固体弹性变形问题的优势，并且在两者耦合过程中进行合理的信息交换，保证计算的快捷性和计算结果的准确性，并且能够自然方便地捕捉流体界面及模拟流体形态，无需花很大精力对网格进行复杂操作来保证计算顺利进行。对网格进行复杂操作来保证计算顺利进行。对网格进行复杂操作来保证计算顺利进行。

【技术实现步骤摘要】
一种适用于模拟水弹性问题的耦合数值计算方法


[0001]本专利技术涉及计算力学领域，尤其是涉及一种适用于模拟水弹性问题的耦合数值计算方法。

技术介绍

[0002]流固耦合现象在自然界和工业生产中都十分常见，例如海鸥水中捕鱼，民航飞机水上迫降、轮船水上航行等等，涉及到砰击、自由液面演化和入水空泡等复杂物理现象，具有强非线性特点，吸引了很多科研工作者的关注。对于流固耦合过程的模拟预测，理论研究具有较大的局限性，难以给出复杂结构作用过程的解析解，通常仅能常用于解决简单外形结构的流固作用。目前在该领域中，主要采用实验研究和数值研究两种方法。其中实验研究是最直观可信的方案，研究人员采用实验方法在固体受力载荷预报和多相流动等方面做了大量的工作并给出了直观可信的结果，但缺点是对试验设施要求较高，且花费成本相对较大。
[0003]随着近些年计算机硬件水平和计算方法的不断飞速发展，数值计算越来越多地被用于模拟流固耦合问题。目前数值计算方法主要分为基于网格的数值方法和基于粒子的数值方法。基于网格的数值方法主要包括有限差分法(FDM)、有限体积法(FVM)和有限元法(FEM)等，这些方法从计算力学学科发展初始便开始发展，至今方法成熟度已经很高，被广泛应用于固体、流体、电磁等各类领域。粒子类数值计算方法发展时间还较短，主要从二十世纪中后期开始发展，目前应用较为广泛的主要有光滑粒子流体动力学方法(SPH)、再生核粒子法(RKPM)和点插值法(PIM)等等。粒子类方法在自由表面流动，运动边界和大变形等问题的模拟上相比基于网格的方法具有较大优势，因此近些年发展很快。
[0004]对于流固耦合的数值模拟无论是采用基于网格的数值计算方法还是基于粒子的数值计算方法，为图方便目前大多会将固体视为刚体，忽略固体的弹性变形，达不到真实模拟的效果。对于少数考虑结构响应的流固耦合模拟，基于网格数值模拟方法在处理流体自由界面及大变形时需要对网格进行很复杂操作规避网格异变甚至破碎，基于粒子的数值方法则在处理结构响应时往往会遇到较大困难。

技术实现思路

[0005]本专利技术着力解决的问题是流固耦合作用过程，尤其是涉及自由界面的流体弹性体相互作用过程。旨在提供一种实施方便、计算准确的用于模拟流场压力、液面演化和弹性体受力、变形的粒子/网格耦合数值计算方法，为模拟自由表面流动和砰击等水弹性问题提供支撑。本专利技术提出SPH
‑
FEM粒子/网格耦合数值方法来实现对流体弹性体间相互作用的准确模拟。其中SPH方法用来模拟流场，发挥粒子类方法在模拟涉及自由交界面及界面变形问题上的优势；FEM方法来模拟弹性体，发挥网格类方法处理固体弹性变形问题的优势，并且在两者耦合过程中进行合理的信息交换，保证计算的快捷性和计算结果的准确性。
[0006]一种适用于模拟水弹性问题的耦合数值计算方法，其特征在于，分别进行SPH计算
和FEM计算，并进行信息交互：
[0007]S1.采用SPH方法获取流场信息；
[0008]S11.将流体离散成一系列粒子，每个粒子代表一个流体微元；
[0009]令t
f
＝t0，其中，t0为初始时刻；
[0010]S12.计算流体粒子密度；
[0011]S13.计算流场压力；
[0012]S14.根据质量守恒方程和动量守恒方程计算密度变化率和加速度；
[0013]S15.根据运动学方程获得粒子位移；
[0014]S16.令t
f
＝t
f
+Δt
f
，返回步骤S12，其中，Δt
f
为第一时间步长；
[0015]S2.采用FEM计算弹性体结构响应；
[0016]S21.将弹性体离散成一系列网格单元，所述网格单元采用退化的网格单元，令t
e
＝t0，其中，t0为初始时刻；
[0017]S22.计算内应力；
[0018]S23.计算外部载荷；
[0019]S24.根据动力学方程计算加速度；
[0020]S25.根据运动学方程获得弹性体网格点位移；
[0021]S26.令t
e
＝t
e
+Δt
e
，返回步骤S22，其中，Δt
e
为第二时间步长；
[0022]S3.SPH、FEM计算信息交互
[0023]每经过一个第一时间步长，SPH和FEM进行一次信息交互，包括：
[0024]将SPH计算得到的流场压力和流体粒子密度带入步骤S23中计算外部载荷，以及S13计算交界面处的流体粒子受力时，接收来自弹性体的运动边界信息。
[0025]进一步地，S13计算交界面处的流体粒子受力时，接收来自弹性体的运动边界信息包括：
[0026]对于特定流体粒子，若其支持域内有弹性体单元，其受力来源为弹性体单元和流体粒子的共同作用力：
[0027][0028]若其支持域完全处于流体内部，只有流体粒子的作用力：
[0029][0030]其中，为梯度算子，n
k
为单元法向向量，N
f
为支持域内流体粒子总数，N
e
为支持域内弹性体单元总数，W为核函数，S
k
为弹性体单元面积，V
j
为流体粒子体积，r泛指流体计算中的参数。
[0031]进一步地，所述第一时间步长Δt
f
和第二时间步长Δt
e
根据如下方式确定：
[0032][0033]h、c0、v
max
和a
max
分别代表流体粒子的光滑长度、人工声速、最大速度和最大加速度。E和ρ分别是弹性体的弹性模量和密度，Δx为弹性体单元的网格尺寸。
[0034]进一步地，Δt
f
＝nΔt
e
，n为正整数。
[0035]进一步地，步骤S13中，流场压力采用弱可压假设，通过状态方程计算：
[0036][0037]其中，p
f
表示流体粒子压力，ρ
f
表示流体粒子当前密度，ρ
f0
表示流体粒子初始密度，B为常数。
[0038]进一步地，步骤S14中，通过下式计算流体粒子的密度变化率和加速度：
[0039][0040]其中，m表示粒子质量，v
ij
表示流体粒子i和流体粒子j的速度差，流体粒子j是流体粒子i支持域内的粒子，α和β表示坐标方向，W表示核函数，表示对流体粒子i支持域内的N个粒子的相关参数进行叠加求和，x表示位移矢量。
[0041]进一步地，步骤S22中，所述内应力P的计算为：
[0042][0043][0044]其中，N
I
表示一维形函数，b是网格单元的高度，α和β表示坐标方向，x表示位移矢量，X表示初始时刻的位移，ζ是主纤维矢量方向的参数坐标，y...

【技术保护点】

【技术特征摘要】
1.一种适用于模拟水弹性问题的耦合数值计算方法，其特征在于，分别进行SPH计算和FEM计算，并进行信息交互：S1.采用SPH方法获取流场信息；S11.将流体离散成一系列粒子，每个粒子代表一个流体微元；令t
f
＝t0，其中，t0为初始时刻；S12.计算流体粒子密度；S13.计算流场压力；S14.根据质量守恒方程和动量守恒方程计算密度变化率和加速度；S15.根据运动学方程获得粒子位移；S16.令t
f
＝t
f
+Δt
f
，返回步骤S12，其中，Δt
f
为第一时间步长；S2.采用FEM计算弹性体结构响应；S21.将弹性体离散成一系列网格单元，所述网格单元采用退化的网格单元，令t
e
＝t0，其中，t0为初始时刻；S22.计算内应力；S23.计算外部载荷；S24.根据动力学方程计算加速度；S25.根据运动学方程获得弹性体网格点位移；S26.令t
e
＝t
e
+Δt
e
，返回步骤S22，其中，Δt
e
为第二时间步长；S3.SPH、FEM计算信息交互每经过一个第一时间步长，SPH和FEM进行一次信息交互，包括：将SPH计算得到的流场压力和流体粒子密度带入步骤S23中计算外部载荷，以及S13计算交界面处的流体粒子受力时，接收来自弹性体的运动边界信息。2.根据权利要求1所述的一种适用于模拟水弹性问题的耦合数值计算方法，其特征在于，S13计算交界面处的流体粒子受力时，接收来自弹性体的运动边界信息包括：对于特定流体粒子，若其支持域内有弹性体单元，其受力来源为弹性体单元和流体粒子的共同作用力：若其支持域完全处于流体内部，只有流体粒子的作用力：其中，为梯度算子，n
k
为单元法向向量，N
f
为支持域内流体粒子总数，N
e
为支持域内弹性体单元总数，W为核函数，S
k
为弹性体单元面积，V
j
为流体粒子体积，r泛指流体计算中的参数。3.根据权利要求1所述的一种适用于模拟水弹性问题的耦合数值计算方法，其特征在于，所述第一时间步长Δt
f
和第二时间步长Δt
e
根据如下方式确定：
h、c0、v
max
和a
max
分别代表流体粒子的光滑长度、人工声速、最大速度和最大加速度。E和ρ分别是弹性体的弹性模量和密度，Δx为弹性体单元的网格尺寸。4.根据权利要求3所述的一种适用于模拟水...

【专利技术属性】
技术研发人员：施文奎，陈程，沈雁鸣，鞠欣洋，陈坚强，郭勇颜，万云博，
申请(专利权)人：中国空气动力研究与发展中心计算空气动力研究所，
类型：发明
国别省市：