CN116992747B - 一种基于sph流固耦合的冲击式水轮机动力学分析方法 - Google Patents
一种基于sph流固耦合的冲击式水轮机动力学分析方法 Download PDFInfo
- Publication number
- CN116992747B CN116992747B CN202311267496.4A CN202311267496A CN116992747B CN 116992747 B CN116992747 B CN 116992747B CN 202311267496 A CN202311267496 A CN 202311267496A CN 116992747 B CN116992747 B CN 116992747B
- Authority
- CN
- China
- Prior art keywords
- sph
- fluid
- particles
- particle
- density
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000010168 coupling process Methods 0.000 title claims abstract description 83
- 230000008878 coupling Effects 0.000 title claims abstract description 81
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 81
- 239000007787 solid Substances 0.000 title claims abstract description 49
- 238000004458 analytical method Methods 0.000 title claims abstract description 28
- 239000002245 particle Substances 0.000 claims abstract description 245
- 239000012530 fluid Substances 0.000 claims abstract description 114
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 113
- 238000004364 calculation method Methods 0.000 claims abstract description 61
- 238000000034 method Methods 0.000 claims abstract description 57
- 239000006185 dispersion Substances 0.000 claims abstract description 21
- 238000004088 simulation Methods 0.000 claims abstract description 21
- 230000001808 coupling effect Effects 0.000 claims abstract description 10
- 230000010354 integration Effects 0.000 claims description 9
- 238000009792 diffusion process Methods 0.000 claims description 8
- 238000003860 storage Methods 0.000 claims description 8
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000005059 solid analysis Methods 0.000 claims description 3
- 239000007788 liquid Substances 0.000 abstract description 11
- 230000001133 acceleration Effects 0.000 abstract description 7
- 238000005516 engineering process Methods 0.000 abstract description 6
- 230000006870 function Effects 0.000 description 21
- 238000010586 diagram Methods 0.000 description 10
- 238000006243 chemical reaction Methods 0.000 description 4
- 238000004590 computer program Methods 0.000 description 4
- 238000013461 design Methods 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000014509 gene expression Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 239000008187 granular material Substances 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000000889 atomisation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于SPH流固耦合的冲击式水轮机动力学分析方法，方法通过水轮机的网格模型确定SPH颗粒离散；确定水轮机的计算域，并根据实际水流入口确定SPH流体速度入口；确定流体与水轮机的若干耦合参数，耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；根据SPH颗粒离散确定初始条件，根据SPH流体速度入口确定入口边界条件，根据各耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与水轮机的耦合作用仿真求解。本发明中SPH颗粒可直接表征自由液面，无需重新划分网格，并且适合采用GPU并行加速技术，因此本发明可以更好地模拟水轮机与流体之间的耦合运动，并大大提升计算效率。
Description
技术领域
本发明涉及水利水电领域，尤其涉及的是一种基于SPH流固耦合的冲击式水轮机动力学分析方法。
背景技术
水轮发电机是以水轮机为原动机将水能转化为电能的发电机，冲击式水轮机为水轮机中重要机型，其特点是使用水头高，结构简单。目前国产冲击式水轮机转轮型号少且效率偏低，特别是1000米以上水头段的冲击式水轮机设计制造仍是空白。因此，对高效率冲击式水轮机的研究迫在眉睫。
在过去的研究中，试验研究为主要手段，然而由于射流对水轮机的高速冲击过程伴随着强烈的雾化现象，因此通常难以观测到水轮机的内部流态。近年来，随着计算机技术的发展，使得对于冲击式水轮机的整机数值模拟成为可能，并日益受到众多研究者的重视。冲击式水轮机运动属于三维瞬态流固耦合问题，水和水斗之间的相互作用错综复杂，其内部流态及自由液面也非常复杂。采用传统的网格方法分析水轮机与流体之间的耦合运动时，面临着自由液面捕捉困难、不断的网格重构造成计算效率低等难题，难以对冲击式水轮机的全工作历程进行直接模拟。
因此，现有技术还有待改进和发展。
发明内容
本发明要解决的技术问题在于，针对现有技术的上述缺陷，提供一种基于SPH流固耦合的冲击式水轮机动力学分析方法，旨在解决现有技术中采用传统的网格方法分析水轮机与流体之间的耦合运动时，自由液面捕捉困难，并且不断的网格重构会造成计算效率低，因此难以对水轮机的全工作历程进行直接模拟的问题。
本发明解决问题所采用的技术方案如下：
第一方面，本发明实施例提供一种基于SPH流固耦合的冲击式水轮机动力学分析方法，其中，所述方法包括：
通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解。
在一种实施方式中，所述通过水轮机的网格模型确定所述水轮机的SPH颗粒离散，包括：
从所述网格模型的表面由外向内生成预设层数的SPH颗粒，得到所述水轮机对应的SPH颗粒离散，其中，所述预设层数大于或者等于三层。
在一种实施方式中，所述水轮机对应的SPH颗粒与流体SPH颗粒的SPH粒径相同，SPH粒径基于所述水轮机的最小特征尺寸确定。
在一种实施方式中，所述确定所述水轮机对应的计算域，包括：
获取所述水轮机的外形数据和预设的射流扩散范围；
根据所述外形数据和所述射流扩散范围，确定所述计算域。
在一种实施方式中，所述时间积分的计算方式包括：
针对首个时间步长，该时间步长结束时，流体SPH颗粒的密度、速度以及内能由初始状态向前推进半个时间步长，流体SPH颗粒的位置向前推进一个时间步长；
针对非首个的每一时间步长，该时间步长开始时，密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值；
当一个时间步长结束时，粒子的密度、速度、内能和位置向前推进一个时间步长。
在一种实施方式中，所述SPH核函数为：
；
其中，q为粒子距离与光滑长度比值，σ为归一化系数。
在一种实施方式中，所述密度积分的计算方法包括：
；
其中，h为光滑长度，c0为人工声速，ri为粒子i坐标，rj为粒子j坐标，m为粒子质量，ρ为粒子密度，vij为粒子相对速度，δ为耗散系数，Ψij为粒子密度差，为密度耗散项。
在一种实施方式中，所述粒子间排斥力的计算方法包括：
；
其中，为粒子i、j之间的排斥力，n为粒子法向量，v为粒子速度，Wij为核函数，c为人工声速，rij为粒子i、j的距离，h为光滑长度。
第二方面，本发明实施例还提供一种基于SPH流固耦合的冲击式水轮机动力学分析装置，所述装置包括：
固体分析模块，用于通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
流体分析模块，用于确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
参数确定模块，用于确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
仿真求解模块，用于根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解。
第三方面，本发明实施例还提供一种计算机可读存储介质，其上存储有多条指令，其中，所述指令适用于由处理器加载并执行，以实现上述任一所述的基于SPH流固耦合的冲击式水轮机动力学分析方法的步骤。
本发明的有益效果：本发明实施例通过水轮机的网格模型确定水轮机的SPH颗粒离散；确定水轮机的计算域，并根据实际水流入口确定SPH流体速度入口；确定流体与水轮机的若干耦合参数，耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；根据SPH颗粒离散确定初始条件，根据SPH流体速度入口确定入口边界条件，根据各耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与水轮机的耦合作用仿真求解。本发明中SPH颗粒可直接表征自由液面，不需要对网格进行重新划分，并且SPH显示计算方法适合采用GPU并行加速技术，因此通过SPH颗粒可以更好地模拟水轮机与流体之间的耦合运动，并大大提升计算效率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的基于SPH流固耦合的冲击式水轮机动力学分析方法的流程示意图。
图2是本发明实施例提供的水轮机网格示意图及射流入口位置示意图。
图3是本发明实施例提供的水轮机SPH颗粒离散示意图。
图4是本发明实施例提供的射流冲击水轮机流场分布示意图一。
图5是本发明实施例提供的射流冲击水轮机流场分布示意图二。
图6是本发明实施例提供的入口边界条件方法示意图。
图7是本发明实施例提供的水轮机扭矩计算结果示意图。
图8是本发明实施例提供的基于SPH流固耦合的冲击式水轮机动力学分析装置的模块示意图。
图9是本发明实施例提供的终端的原理框图。
具体实施方式
本发明公开了一种基于SPH流固耦合的冲击式水轮机动力学分析方法，为使本发明的目的、技术方案及效果更加清楚、明确，以下参照附图并举实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。 应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或无线耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的全部或任一单元和全部组合。
本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)，具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语，应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样被特定定义，否则不会用理想化或过于正式的含义来解释。
针对现有技术的上述缺陷，本发明提供一种基于SPH流固耦合的冲击式水轮机动力学分析方法，如图1所示，所述方法包括：
步骤S100、通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
步骤S200、确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
步骤S300、确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
步骤S400、根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解。
SPH数值模拟将所有的物质都看作是粒子组成，有利于自由液面捕捉，并且具有无网格特性，无需进行网格重构，可以有效提升计算效率。本实施例的目标是采用SPH颗粒来模拟高速射流冲击水轮机的流固耦合过程，以提高模拟结果的准确性和模拟效率。具体地，如图2所示，本实施例中的网格模型是在水轮机的几何模型表面采用三角面片形成的。网格模型的存储文件格式可以采用stl格式，obj格式，stp格式，step格式中任一一种格式。针对水轮机部分，本实施例会将网格模型离散为合适粒径的拉格朗日粒子，从而得到水轮机的SPH颗粒离散。针对射流部分，本实施例会根据实际水流入口确定SPH流体速度入口。此外本实施例还会预先定义流体与水轮机之间的多个耦合参数，包括但不限于时间积分、SPH核函数、密度积分以及流固耦合。其中，时间积分采用类似蛙跳算法的方法，即流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；SPH核函数则是采用预设样条函数，例如Wendland五次样条函数；密度积分中增加了密度耗散计算，以解决常规SPH方法存在压力场计算不准确的问题；流固耦合采用粒子间排斥力的计算方法，例如Han&Zhang排斥力方法。最后将前述步骤得到的SPH颗粒离散作为初始条件，SPH流体速度入口作为入口边界条件，各耦合参数作为求解参数，采用显式中心差分时间推进的方法，实现流体与水轮机的耦合作用仿真求解。本发明中SPH颗粒可直接表征自由液面，不需要对网格进行重新划分，并且SPH显示计算方法适合采用GPU并行加速技术，因此本发明通过SPH颗粒可以更好地模拟水轮机与流体之间的耦合运动，并大大提升计算效率。
在一种实现方式中，所述通过水轮机的网格模型确定所述水轮机的SPH颗粒离散，包括：
从所述网格模型的表面由外向内生成预设层数的SPH颗粒，得到所述水轮机对应的SPH颗粒离散，其中，所述预设层数大于或者等于三层。
具体地，如图3所示，从网格模型的表面出发向内生成预设层数的SPH颗粒，这些离散的SPH颗粒即可反映水轮机的固体形态。为了保证计算精度，水轮机对应的SPH颗粒的层数至少需要3层。
在一种实现方式中，所述水轮机对应的SPH颗粒与流体SPH颗粒的SPH粒径相同，SPH粒径基于所述水轮机的最小特征尺寸确定。
具体地，水轮机对应的SPH颗粒与流体SPH颗粒的SPH粒径相同，均取决于水轮机的网格模型的最小特征尺寸。以水轮机水斗的厚度为例，为保证沿水斗厚度方向有三层SPH颗粒，SPH粒径应为厚度的1/3。
在一种实现方式中，所述根据实际水流入口确定SPH流体速度入口，包括：
根据实际水流入口确定若干参数，得到所述SPH流体速度入口，其中，所述参数包括入口流量、直径以及流体密度。
具体地，水轮机在实际射流位置可以反映实际水流入口（如图2所示），本实施例需要在实际水流入口处设置SPH入口边界条件，定义射流位置、入口流量、直径、流体密度等参数，从而得到SPH流体速度入口。例如，射流直径为156mm，射流速度为78.5m/s，流体密度为1000kg/m3。
在一种实现方式中，所述确定所述水轮机对应的计算域，包括：
获取所述水轮机的外形数据和预设的射流扩散范围；
根据所述外形数据和所述射流扩散范围，确定所述计算域。
具体地，本实施例主要是研究计算域内的流固耦合运动，所以计算域至少应包含整个水轮机。并且为保证射流冲击水斗后流场的准确计算，还需要预估射流扩散范围，并使得计算域包含射流扩散范围。图4、5展示了射流冲击水轮机流场分布示意图。
在一种实现方式中，所述时间积分的计算方式包括：
针对首个时间步长，该时间步长结束时，流体SPH颗粒的密度、速度以及内能由初始状态向前推进半个时间步长，流体SPH颗粒的位置向前推进一个时间步长；
针对非首个的每一时间步长，该时间步长开始时，密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值；
当一个时间步长结束时，粒子的密度、速度、内能和位置向前推进一个时间步长。
具体地，本实施例中采用蛙跳算法进行时间积分，在第一个时间步长结束时，密度、速度和内能由初始状态（可以由用户设定的流体初始速度和密度确定）向前推进半个时间步长，而粒子的位置向前推进一个时间步长：
；
；
；
其中，ρ为密度，t为时间，v为速度，x为位移；
在之后的每一个时间步的开始，粒子的密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值，以便与粒子的位置时刻保持一致：
；
；
当一个时间步长结束时，粒子的密度、速度、内能和位置向前推进一个时间步长：
；
；
。
在一种实现方式中，所述SPH核函数为：
；
其中，q为粒子距离与光滑长度比值，σ为归一化系数，在一维，二维，三维问题中分别可以取2/3，10/（7π）和1/π。
在一种实现方式中，所述密度积分的计算方法包括：
；
其中，h为光滑长度，c0为人工声速，ri为粒子i的坐标，rj为粒子j的坐标，m为粒子质量，ρ为粒子密度，vij为粒子相对速度，δ为耗散系数，Ψij为粒子密度差，为密度耗散项。δ取值可以为0.1，Ψij=ρ j -ρ i 。
具体地，常规的SPH方法计算得到的密度含有高频噪点，导致流体的压力场容易产生紊乱。因此为了解决压力场计算不准确的问题，本实施例提供了一种类似半δ-SPH方法，在常规SPH方法的离散连续性方程中添加了一个密度耗散项。
在一种实现方式中，粒子密度差的计算方法为：
；
其中，为粒子i的修正密度梯度。本实施例在计算粒子密度差时加入了修正密度梯度，可以得到更加准确的结果。
在一种实现方式中，所述粒子间排斥力的计算方法包括：
；
其中，为粒子i、j之间的排斥力，n为粒子法向量，v为粒子速度，W为核函数，c为人工声速，rij为粒子i、j的距离，h为光滑长度。
在一种实现方式中，基于SPH流固耦合的冲击式水轮机的动量方程求解方法，包括：
当采用人工粘性时，离散的动量方程为：
；
其中，右边第一项为压力引起的加速度，压力的计算采用状态方程，其表达式为：，/>为参考密度。第二项为粘性耗散项，/>为人工粘性系数，一般取为0.1，的表达式为：/>。
人工粘性方式的动量方程可以用于处理无粘流问题，此时的人工粘性是能量的耗散项。
在另一种实现方式中，当采用人工粘性与真实粘性时，首先确定人工粘性与真实粘性之间的对应关系：
；
其中，d是求解问题的维度。
在另一种实现方式中，针对低雷诺数管道流采用的加速度计算公式为：
；
其中，为i粒子的真实粘性系数，/>为粒子间距离。
在一种实现方式中，初始条件包括粒子的速度场与压力场的初始赋值。针对已知的初始速度和压力分布，通过前端输入表达式（例如v=v（x，y，z））即可完成初始条件的设定。
在一种实现方式中，流体与水轮机之间采用无滑移边界条件。
具体地，可以设定固壁粒子的速度为内部流体粒子外插后直接反向，或者边界粒子平行于壁面的速度与内部流体粒子速度外插值相反，垂直于壁面的速度与内部流体粒子速度外插值也相反，以实现无滑移边界条件：
；
；
其中，为自定义墙体运动，/>为流体粒子插值速度，/>为边界虚粒子。边界粒子压力插值如下：
。
在一种实现方式中，入口边界条件方法，包括：
在入口区域设置第一缓冲区，其中，第一缓冲区内的流体SPH颗粒根据预设的流速动态生成，第一缓冲区由域边缘和缓冲阈值组成；
将第一缓冲区内的流体SPH颗粒镜像至流体区域内，得到若干镜像粒子；
根据1阶精度插值方法，将镜像粒子的压力值外插至第一缓冲区，得到第一缓冲区内的流体SPH颗粒的压力。
具体地，如图6所示，在入口区域设置第一缓冲区（虚线之间），第一缓冲区的宽度一般大于或等于支持域的长度。该区域的粒子根据给定流速动态生成，粒子速度一般为给定值，但是粒子的压力未知。由于内部流体粒子的压力是已知的，所以本实施例将第一缓冲区粒子镜像到流体区域，得到镜像粒子。然后使用1阶精度插值方法，将镜像粒子的压力值外插至第一缓冲区，从而完成第一缓冲区粒子的压力计算。
举例说明，假设第一缓冲区粒子压力为p o ，位置矢量为r 0 ，镜像粒子的压力为p k ，位置矢量为r k ，镜像粒子处的修正核函数梯度为，则第一缓冲区粒子的压力计算公式为：
。
在一种实现方式中，出口边界条件方法包括：
在出口区域设置第二缓冲区，当流体SPH颗粒流入第二缓冲区后，平行于流动方向的速度不变，而垂直于流动方向的流速变为0，压力维持不变。
具体地，在出口区域，如果进行简单的粒子删除，将产生强烈的边界反射效应，影响内部流体的计算精度。因此需要设置第二缓冲区，第二缓冲区的宽度一般大于或等于支持域的长度。当粒子流入缓冲区后，平行于流动方向的速度不变，而垂直于流动方向的流速变为0，压力也维持不变，这样在出口区域，其速度梯度和压力梯度为0，从而削弱了边界的反射效应。
本实施例提供的入口边界方法和出口边界方法可以减少出入口边界反射现象，解决出入口不稳定性问题，提高了计算精度。
在一种实现方式中，所述方法还包括：
对所述水轮机的每一SPH颗粒的扭矩求和，得到所述水轮机的整体扭矩。
具体地，如图7所示，本实施例还可以基于每一SPH颗粒的扭矩计算出水轮机的整体扭矩，进而为水轮机的优化设计提供参考数据。
在一种实现方式中，所述方法还包括：
根据所述水轮机的扭矩、转速以及圆周率，确定所述水轮机对应的功率。
具体地，本实施例还可以基于水轮机的扭矩和转速，计算出水轮机的功率，计算公式如下所示：
；
其中，P为水轮机功率，为圆周率，T为水轮机扭矩，n为水轮机转速。
在一种实现方式中，所述方法还包括：
根据所述水轮机的功率、时间、射流质量以及射流速度，确定所述水轮机的能量转换效率。
具体地，本实施例还可以计算水轮机的能量转换效率，通过模拟流固耦合过程实现水轮机动力学分析。计算公式如下所示：
；
其中，为能量转换效率，t为时间，m为射流质量，v为射流速度。
本发明的优点在于：
1.冲击式水轮机运作时水和水斗之间的相互作用错综复杂。传统的网格方法处理这类问题时，面临着自由液面捕捉困难、不断的网格重构造成的计算效率低等难题，难以对冲击式水轮机的全工作历程进行直接模拟。
2.本发明采用SPH颗粒对流体和水轮机进行离散，SPH颗粒可直接携带复杂流场信息，有效解决了传统网格方法自由液面捕捉困难的问题。
3.本发明可以更好地处理自由液面大变形问题，便于水轮机与流体之间的耦合运动计算，得到水轮机内部流态，并能直接提取水轮机扭矩及能量转换效率。
4.SPH计算适合采用GPU并行加速技术，从而便于针对实际工程问题开展大规模三维瞬态数值模拟，为冲击式水轮机的优化设计提供了新的研究思路。
基于上述实施例，本发明还提供了一种基于SPH流固耦合的冲击式水轮机动力学分析装置，如图8所示，所述装置包括：
固体分析模块01，用于通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
流体分析模块02，用于确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
参数确定模块03，用于确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
仿真求解模块04，用于根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解。
基于上述实施例，本发明还提供了一种终端，其原理框图可以如图9所示。该终端包括通过系统总线连接的处理器、存储器、网络接口、显示屏。其中，该终端的处理器用于提供计算和控制能力。该终端的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该终端的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现基于SPH流固耦合的冲击式水轮机动力学分析方法。该终端的显示屏可以是液晶显示屏或者电子墨水显示屏。
本领域技术人员可以理解，图9中示出的原理框图，仅仅是与本发明方案相关的部分结构的框图，并不构成对本发明方案所应用于其上的终端的限定，具体的终端可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一种实现方式中，所述终端的存储器中存储有一个以上的程序，且经配置以由一个以上处理器执行所述一个以上程序包含用于进行基于SPH流固耦合的冲击式水轮机动力学分析方法的指令。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本发明所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。
综上所述，本发明公开了一种基于SPH流固耦合的冲击式水轮机动力学分析方法，所述方法通过水轮机的网格模型确定水轮机的SPH颗粒离散；确定水轮机的计算域，并根据实际水流入口确定SPH流体速度入口；确定流体与水轮机的若干耦合参数，耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；根据SPH颗粒离散确定初始条件，根据SPH流体速度入口确定入口边界条件，根据各耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与水轮机的耦合作用仿真求解。本发明中SPH颗粒可直接表征自由液面，不需要对网格进行重新划分，并且SPH显示计算方法适合采用GPU并行加速技术，因此通过SPH颗粒可以更好地模拟水轮机与流体之间的耦合运动，并大大提升计算效率。
应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (9)
1.一种基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述方法包括：
通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解；
所述时间积分的计算方式为：在第一个时间步长结束时，流体SPH颗粒的密度、速度和内能由初始状态向前推进半个时间步长，而流体SPH颗粒的位置向前推进一个时间步长：
；
；
；
其中，ρ为粒子密度，t为时间，v为速度，x为位移，为第i个粒子的粒子密度，/>为第i个粒子在1/2时间步的粒子密度，/>为时间步长，/>为第i个粒子在初始时间步的粒子密度，/>为第i个粒子，/>为第i个粒子的速度，/>为第i个粒子在初始时间步的速度，/>为第i个粒子在1/2时间步的速度，/>为第i个粒子在第1时间步的位置，/>为第i个粒子在初始时间步的位置；
在之后的每一个时间步的开始，流体SPH颗粒的密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值：
；
；
其中，为第i个粒子在第n个时间步的粒子密度，/>为第i个粒子在第n个时间步的速度；
当一个时间步长结束时，流体SPH颗粒的密度、速度、内能和位置向前推进一个时间步长：
；
；
；
所述密度积分的计算方法包括：
；
其中，h为光滑长度，c0为人工声速，ri为粒子i的坐标，rj为粒子j的坐标，m为粒子质量，ρ为粒子密度，vij为粒子相对速度，δ为耗散系数，Ψij为粒子密度差，为密度耗散项，/>为第j个粒子的粒子质量，/>为第j个粒子的粒子密度，/>为第i、j个粒子之间的核函数梯度；
所述粒子密度差的计算方法为：
；
其中，为粒子i的修正密度梯度，/>为粒子j的修正密度梯度；
所述方法还包括：
在入口区域设置第一缓冲区，其中，所述第一缓冲区内的流体SPH颗粒根据预设的流速动态生成，所述第一缓冲区由域边缘和缓冲阈值组成；
将所述第一缓冲区内的流体SPH颗粒镜像至流体区域内，得到若干镜像粒子；
根据1阶精度插值方法，将所述镜像粒子的压力值外插至所述第一缓冲区，得到所述第一缓冲区内的流体SPH颗粒的压力；
在出口区域设置第二缓冲区，当流体SPH颗粒流入第二缓冲区后，平行于流动方向的速度不变，而垂直于流动方向的流速变为0，压力维持不变。
2.根据权利要求1所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述通过水轮机的网格模型确定所述水轮机的SPH颗粒离散，包括：
从所述网格模型的表面由外向内生成预设层数的SPH颗粒，得到所述水轮机对应的SPH颗粒离散，其中，所述预设层数大于或者等于三层。
3.根据权利要求2所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述水轮机对应的SPH颗粒与流体SPH颗粒的SPH粒径相同，SPH粒径基于所述水轮机的最小特征尺寸确定。
4.根据权利要求1所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述确定所述水轮机对应的计算域，包括：
获取所述水轮机的外形数据和预设的射流扩散范围；
根据所述外形数据和所述射流扩散范围，确定所述计算域。
5.根据权利要求1所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述时间积分的计算方式包括：
针对首个时间步长，该时间步长结束时，流体SPH颗粒的密度、速度以及内能由初始状态向前推进半个时间步长，流体SPH颗粒的位置向前推进一个时间步长；
针对非首个的每一时间步长，该时间步长开始时，密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值；
当一个时间步长结束时，粒子的密度、速度、内能和位置向前推进一个时间步长。
6.根据权利要求1所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述SPH核函数为：
；
其中，q为粒子距离与光滑长度比值，σ为归一化系数。
7.根据权利要求1所述的基于SPH流固耦合的冲击式水轮机动力学分析方法，其特征在于，所述粒子间排斥力的计算方法包括：
；
其中，为粒子i、j之间的排斥力，n为粒子法向量，v为粒子速度，W为核函数，c为人工声速，rij为粒子i、j的距离，h为光滑长度，/>为第i个粒子的粒子速度，/>为第j个粒子的粒子速度，/>为第j个粒子的粒子法向量，/>为核函数。
8.一种基于SPH流固耦合的冲击式水轮机动力学分析装置，其特征在于，所述装置包括：
固体分析模块，用于通过水轮机的网格模型确定所述水轮机的SPH颗粒离散；
流体分析模块，用于确定所述水轮机对应的计算域，并根据实际水流入口确定SPH流体速度入口；
参数确定模块，用于确定流体与所述水轮机的若干耦合参数，其中，所述耦合参数包括：时间积分、SPH核函数、密度积分以及流固耦合；所述时间积分中流体SPH颗粒的密度、速度、内能与流体SPH颗粒的位置分别采用不同的推进方式；所述SPH核函数采用预设样条函数；所述密度积分中包括密度耗散计算；所述流固耦合采用粒子间排斥力的计算方法；
仿真求解模块，用于根据所述SPH颗粒离散确定初始条件，根据所述SPH流体速度入口确定入口边界条件，根据各所述耦合参数确定求解参数；通过显式中心差分时间推进的方法进行流体与所述水轮机的耦合作用仿真求解；
所述时间积分的计算方式为：在第一个时间步长结束时，流体SPH颗粒的密度、速度和内能由初始状态向前推进半个时间步长，而流体SPH颗粒的位置向前推进一个时间步长：
；
；
；
其中，ρ为粒子密度，t为时间，v为速度，x为位移，为第i个粒子的粒子密度，/>为第i个粒子在1/2时间步的粒子密度，/>为时间步长，/>为第i个粒子在初始时间步的粒子密度，/>为第i个粒子，/>为第i个粒子的速度，/>为第i个粒子在初始时间步的速度，/>为第i个粒子在1/2时间步的速度，/>为第i个粒子在第1时间步的位置，/>为第i个粒子在初始时间步的位置；
在之后的每一个时间步的开始，流体SPH颗粒的密度、速度和内能向前再推进半个时间步长，获得整数时间步上的值：
；
；
其中，为第i个粒子在第n个时间步的粒子密度，/>为第i个粒子在第n个时间步的速度；
当一个时间步长结束时，流体SPH颗粒的密度、速度、内能和位置向前推进一个时间步长：
；
；
；
所述密度积分的计算方法包括：
；
其中，h为光滑长度，c0为人工声速，ri为粒子i的坐标，rj为粒子j的坐标，m为粒子质量，ρ为粒子密度，vij为粒子相对速度，δ为耗散系数，Ψij为粒子密度差，为密度耗散项，/>为第j个粒子的粒子质量，/>为第j个粒子的粒子密度，/>为第i、j个粒子之间的核函数梯度；
所述粒子密度差的计算方法为：
；
其中，为粒子i的修正密度梯度，/>为粒子j的修正密度梯度；
所述方法还包括：
在入口区域设置第一缓冲区，其中，所述第一缓冲区内的流体SPH颗粒根据预设的流速动态生成，所述第一缓冲区由域边缘和缓冲阈值组成；
将所述第一缓冲区内的流体SPH颗粒镜像至流体区域内，得到若干镜像粒子；
根据1阶精度插值方法，将所述镜像粒子的压力值外插至所述第一缓冲区，得到所述第一缓冲区内的流体SPH颗粒的压力；
在出口区域设置第二缓冲区，当流体SPH颗粒流入第二缓冲区后，平行于流动方向的速度不变，而垂直于流动方向的流速变为0，压力维持不变。
9.一种计算机可读存储介质，其上存储有多条指令，其特征在于，所述指令适用于由处理器加载并执行，以实现上述权利要求1-7任一所述的基于SPH流固耦合的冲击式水轮机动力学分析方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311267496.4A CN116992747B (zh) | 2023-09-28 | 2023-09-28 | 一种基于sph流固耦合的冲击式水轮机动力学分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311267496.4A CN116992747B (zh) | 2023-09-28 | 2023-09-28 | 一种基于sph流固耦合的冲击式水轮机动力学分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116992747A CN116992747A (zh) | 2023-11-03 |
| CN116992747B true CN116992747B (zh) | 2024-03-22 |
Family
ID=88530681
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311267496.4A Active CN116992747B (zh) | 2023-09-28 | 2023-09-28 | 一种基于sph流固耦合的冲击式水轮机动力学分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116992747B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117252131B (zh) * | 2023-11-20 | 2024-03-01 | 深圳十沣科技有限公司 | 一种适用于薄壁结构物的数值仿真方法及装置 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112711882A (zh) * | 2020-12-29 | 2021-04-27 | 重庆建筑科技职业学院 | 一种冲击式水轮机转轮疲劳破坏模型构建方法 |
| CN113361156A (zh) * | 2021-05-20 | 2021-09-07 | 南京航空航天大学 | 一种基于sph方法的飞机水箱投汲水数值模拟方法 |
| CN113947003A (zh) * | 2021-10-15 | 2022-01-18 | 西安交通大学 | 一种面向热流耦合场景的粒子型无网格仿真系统 |
| CN114662425A (zh) * | 2022-05-25 | 2022-06-24 | 浙江远算科技有限公司 | 一种水轮机启停工况流场仿真预测方法及系统 |
| CN115544663A (zh) * | 2022-10-17 | 2022-12-30 | 南京航空航天大学 | 一种考虑柔性弹体多舱段连接的弹体入水数值模拟方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5670832B2 (ja) * | 2011-05-24 | 2015-02-18 | 富士通株式会社 | シミュレーション方法、シミュレーション装置及びシミュレーションプログラム |
-
2023
- 2023-09-28 CN CN202311267496.4A patent/CN116992747B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112711882A (zh) * | 2020-12-29 | 2021-04-27 | 重庆建筑科技职业学院 | 一种冲击式水轮机转轮疲劳破坏模型构建方法 |
| CN113361156A (zh) * | 2021-05-20 | 2021-09-07 | 南京航空航天大学 | 一种基于sph方法的飞机水箱投汲水数值模拟方法 |
| CN113947003A (zh) * | 2021-10-15 | 2022-01-18 | 西安交通大学 | 一种面向热流耦合场景的粒子型无网格仿真系统 |
| CN114662425A (zh) * | 2022-05-25 | 2022-06-24 | 浙江远算科技有限公司 | 一种水轮机启停工况流场仿真预测方法及系统 |
| CN115544663A (zh) * | 2022-10-17 | 2022-12-30 | 南京航空航天大学 | 一种考虑柔性弹体多舱段连接的弹体入水数值模拟方法 |
Non-Patent Citations (3)
| Title |
|---|
| 冲击式水轮机流固耦合数值模拟研究;蒋勇其;《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》(第8期);第18-45页 * |
| 基于流固耦合的混流式水轮机转轮强度分析;阚阚;郑源;赵连辉;乔木;张策;尉青连;;《水电站机电技术》(第2期);全文 * |
| 航行体水下发射流固耦合效应分析;杜特专;王一伟;黄晨光;廖丽涓;;力学学报(04);第47-57页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116992747A (zh) | 2023-11-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Coutier-Delgosha et al. | Experimental and numerical studies in a centrifugal pump with two-dimensional curved blades in cavitating condition | |
| CN116992747B (zh) | 一种基于sph流固耦合的冲击式水轮机动力学分析方法 | |
| Kang et al. | Experimental study on the three-dimensional flow within a compressor cascade with tip clearance: Part I—Velocity and pressure fields | |
| Pini et al. | Adjoint method for shape optimization in real-gas flow applications | |
| Yang et al. | Toward excellence in turbomachinery computational fluid dynamics: A hybrid structured-unstructured reynolds-averaged navier-stokes solver | |
| CN113111608B (zh) | 一种新型局部湍流脉动生成方法 | |
| CN114357571B (zh) | 建成建筑环境下大气边界层风场特性的反演方法和系统 | |
| CN115906718B (zh) | 一种旋转机械cfd系统 | |
| Shao et al. | Experimental and numerical study of external performance and internal flow of a molten salt pump that transports fluids with different viscosities | |
| Immonen | 2D shape optimization under proximity constraints by CFD and response surface methodology | |
| Zhao et al. | Combined experimental and numerical analysis of cavitating flow characteristics in an axial flow waterjet pump | |
| Ye et al. | Assessment of turbulence models for the boundary layer transition flow simulation around a hydrofoil | |
| Li et al. | Research on time sequence prediction of the flow field structure of supersonic cascade channels in wide range based on artificial neural network | |
| Templalexis et al. | Development of a two-dimensional streamline curvature code | |
| Weatherill et al. | Grid adaptation using a distribution of sources applied to inviscid compressible flow simulations | |
| CN113128096B (zh) | 一种获取水下航行器直航附加质量方法 | |
| CN110705189A (zh) | 一种沉淀气浮池气浮区流体力学模型的建立方法 | |
| Liu et al. | Lift-drag characteristics of S-shaped hydrofoil under different cloud cavitation conditions | |
| Kadhim et al. | Design optimization workflow and performance analysis for contoured endwalls of axial turbines | |
| Bao et al. | Smoothed particle hydrodynamics with κ-ε closure for simulating wall-bounded turbulent flows at medium and high Reynolds numbers | |
| Meng et al. | Fast flow prediction of airfoil dynamic stall based on Fourier neural operator | |
| Mannion et al. | A CFD investigation of a variable-pitch vertical axis hydrokinetic turbine with incorporated flow acceleration | |
| Zhang et al. | Research and development on the hydraulic design system of the guide vanes of multistage centrifugal pumps | |
| Wang et al. | Investigation of cavitation and flow characteristics of tip clearance of bidirectional axial flow pump with different clearances | |
| Markus et al. | A Virtual Free Surface (VFS) model for efficient wave–current CFD simulation of fully submerged structures |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |