CN107423511A - 满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 - Google Patents
满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 Download PDFInfo
- Publication number
- CN107423511A CN107423511A CN201710628403.4A CN201710628403A CN107423511A CN 107423511 A CN107423511 A CN 107423511A CN 201710628403 A CN201710628403 A CN 201710628403A CN 107423511 A CN107423511 A CN 107423511A
- Authority
- CN
- China
- Prior art keywords
- boundary
- force
- pressure
- condition
- time
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 60
- 238000004364 calculation method Methods 0.000 claims abstract description 29
- 238000012937 correction Methods 0.000 claims abstract description 8
- 238000007654 immersion Methods 0.000 claims description 26
- 239000006185 dispersion Substances 0.000 claims description 18
- 239000012530 fluid Substances 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 6
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 230000008878 coupling Effects 0.000 abstract description 6
- 238000010168 coupling process Methods 0.000 abstract description 6
- 238000005859 coupling reaction Methods 0.000 abstract description 6
- 238000010586 diagram Methods 0.000 description 8
- 238000007598 dipping method Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000017531 blood circulation Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 210000003709 heart valve Anatomy 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/06—Power analysis or power optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，主要包括以下步骤：采用粘性不可压流N‑S控制方程，将动量方程沿特征线展开进行时间离散，利用有限元法进行空间离散，在离散后的控制方程中加入附加体积力；采用直接力法求解附加体积力，使边界上的插值速度等于期望的速度分布；采用分步法时，将附加体积力按是否与压力耦合分为两部分，分别求解；在附加体积力与压力的迭代计算内进行速度校正，得到最终的流场分布。本发明解决了现有的浸入边界法中，无滑边界条件和连续性条件不能同时满足的问题；将附加体积力分为两部分后，与压力不耦合的分力项的求解成本将有效降低，提高了浸入边界法的计算效率。
Description
技术领域
本发明涉及一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，属于浸入边界法数值求解技术领域。
背景技术
1972年Peskin首先提出浸入边界法，并成功将其应用于生物流体力学，模拟人体心脏瓣膜中的血液流动。其基本思想是将结构的边界模化成Navier-Stokes动量方程中的附加体积力。采用结构化的笛卡尔网格，同时采用两种坐标系，欧拉坐标系描述流场变量，拉格朗日坐标系描述浸入边界变量。
当流场中存在形状复杂的边界时，若采用贴体网格，生成网格的成本较高。在浸入边界法中采用笛卡尔网格，网格的生成更容易。尤其是处理移动边界问题时，对于贴体网格，不仅需要在每一个时间步更新网格，而且需要交换新旧网格中的数据。若采用笛卡尔网格，则不需要更新网格，大大提高了处理动边界问题的效率。
起初，Peskin教授提出的方法主要用于解决弹性边界问题，后来逐步发展为可以解决刚性边界问题。目前，浸入边界法的形式主要分为两种。一是由Lai和Peskin教授提出的反馈力法，该方法通过在动量方程中加入反馈力使得边界上的相对速度为零，相当于在网格交界点上施加了虚拟弹簧与阻尼。但存在很多问题，例如存在较大的虚拟振荡，需要严格限制时间步长等。二是由Mohd-Yusof教授提出的离散直接力法，该方法通过在动量方程中加入附加体积力，直接使流场插值速度等于边界上期望的速度分布。该方法不需要采用小时间步，所以节约了计算成本且能够有效解决三维问题。目前直接力法是最好的解法且应用广泛。
在直接力法中，最大的难点在于求解离散N-S方程。最常见的解法是显式方法，即在当前时间步的计算中，直接使用上一个时间步的压力，没有考虑附加体积力与压力的耦合。另一种解法是隐式方法，即使用当前时间步的压力，计算当前时间步的附加体积力，该方法考虑了附加体积力与压力的耦合。但目前直接力法的算法中，无滑边界条件和连续性条件不能同时满足。及春宁教授提出了一种基于嵌入式迭代的高精度浸入边界法，在压力泊松方程的迭代计算中，嵌入附加体积力的计算。虽然有效提高了施加无滑边界条件的精度，但在附加体积力与压力的迭代计算完成后进行速度校正，将使得结果不能满足无滑边界条件。
发明内容
为解决现有技术的不足，本发明的目的在于提供一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，解决了现有技术中求解离散流体N-S方程时，没有考虑附加体积力与压力的耦合、无滑边界条件和连续性条件不能同时满足的问题。
为了实现上述目标，本发明采用如下的技术方案：
一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，包括如下步骤：
步骤1)采用粘性不可压流N-S控制方程，将动量方程沿特征线展开进行时间离散，利用有限元法进行空间离散，在离散后的控制方程中加入附加体积力；
步骤2)采用直接力法，使得流场插值速度等于边界上期望的速度分布，得到求解式；
步骤3)根据步骤2)的求解式，将附加体积力按是否与压力耦合分为两部分，记为fa、fb；
步骤4)采用分步法求解格式，根据无滑边界条件，对与压力不耦合的附加体积力fa进行迭代计算，得到第n+1个时间步下的faΔt；
步骤5)根据无滑边界条件，对与压力耦合的附加体积力fb进行迭代求解，结合第k个迭代步的fb,(k)Δt和Hn+1,(k)，得到第k+1个迭代步的fb,(k+1)Δt；
步骤6)结合步骤5)中第k+1个迭代步中的fb,(k+1)Δt，根据连续性条件，得到压力第k+1个迭代步的值Hn+1,(k+1)；
步骤7)根据步骤4)-6)中的计算结果，得到第n+1个时间步下、第k+1个迭代步的校正速度un+1,(k+1)；
步骤8)重复步骤4)-7)，迭代计算直至收敛，得到第n+1个时间步的流速un+1；
步骤9)重复步骤4)-8)，进行下一个时间步的迭代计算，直至计算完所有时间步。
进一步地，所述步骤1)具体步骤如下：
步骤1-1)采用粘性不可压流N-S控制方程，具体为：
动量方程：
连续性方程：
其中，为梯度算子，τ表示粘性应力张量，μ表示动力粘滞系数，u表示流体速度，t表示时间，ρ是流体密度，p是压力，fb为重力等体积力。
步骤1-2)将动量方程沿特征线展开，得到时间离散后的动量方程，加入附加体积力，具体为：
其中，Δu为第n+1个和第n个时间步速度的差，Δt是计算时间步长，pn+1为第n+1个时间步的压力，fn+1为第n+1个时间步的附加体积力。
进一步地，引入两种算子Gn和Hn+1将离散后的动量方程重新表示，
其中，Gn为时间离散后的动量方程中，第n个时间步中，与速度un相关的项；Hn+1为时间离散后的动量方程中第n+1个时间步中，与压力pn+1相关的项。
离散后的动量方程可表示为：Δu＝Gn+Hn+1+fn+1Δt，即un+1＝un+Gn+Hn+1+fn+1Δt
其中，un+1为第n+1个时间步的速度。
进一步地，所述步骤2)具体步骤如下：
步骤2-1)定义大写字母代表浸入边界点上的物理量，定义小写字母代表笛卡尔网格交界点上的物理量；
步骤2-2)定义两个函数， 其中，表示在网格交界点上的物理量插值到浸入边界点上，D(Φ)表示在浸入边界点上的物理量投射到网格交界点上，表示在网格交界点上的物理量，Φ表示在浸入边界点上的物理量，u(xi)为网格交界点xi上的速度，U(Xj)为浸入边界点Xj上的速度，δ(xi-Xj)为网格交界点xi与浸入边界点Xj上的物理量进行插值计算时所采用的插值函数，Ω为插值计算中涉及到的网格单元，ΔΩi是网格交界点xi的体积，F(Xj)为浸入边界点Xj上的力，f(xi)为xi网格交界点上的力，ΔVj是浸入边界点Xj的体积，N为浸入边界点个数；
步骤2-3)采用直接力法，使得流场插值速度等于边界上期望的速度分布，即满足无滑边界条件，浸入边界点上的速度等于笛卡尔网格交界点的插值速度：I(un+1)＝Ud，其中，Ud是第n+1个时间步边界上期望的速度分布；
根据动量离散方程，取算子：
I(un+Gn+Hn+1+fn+1Δt)＝Ud，即I(fn+1Δt)＝Ud-I(un+Gn+Hn+1)
D[I(fn+1Δt)]＝D[Ud-I(un+Gn+Hn+1)]
若D[I(fn+1Δt)]＝fn+1Δt，得到求解式fn+1Δt＝D[Ud-I(un+Gn+Hn+1)]。
进一步地，所述步骤3)中将附加体积力fn+1按是否与压力耦合分为两部分：fn+1Δt＝faΔt+fbΔt，其中，fa为第n+1个时间步中与压力不耦合的附加体积力，fb为第n+1个时间步中与压力耦合的附加体积力，faΔt＝D[Ud-I(un+Gn)]，fbΔt＝-D[I(Hn+1)]。
进一步地，所述步骤4)的具体步骤为：
步骤4-1)令fn,(k+1)Δt＝fn,(k)Δt+Δfn,(k)Δt，其中，fn,(k+1)为第n个时间步下，第k+1个迭代步的力，fn,(k)为第n个时间步下，第k个迭代步的力，Δfn,(k)为第n个时间步下，第k个和第k+1个迭代步力的迭代误差；
步骤4-2)根据无滑边界条件：Δfa，(k)Δt＝D[Ud-I(un+Gn+fa，(k)Δt)]，fa,(k+1)Δt＝fa,(k)Δt+Δfa,(k)Δt，其中，Δfa,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力不耦合的附加体积力的迭代误差；fa,(k)为第n+1个时间步下，第k个迭代步，与压力不耦合的附加体积力；fa,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力不耦合的附加体积力；
迭代计算，直至收敛，得到第n+1个时间步faΔt的值。
进一步地，所述步骤5)具体过程如下：
Δfb，(k)Δt＝D[Ud-I(un+Gn+faΔt+fb，(k)Δt+Hn+1，(k))]，fb,(k+1)Δt＝fb,(k)Δt+Δfb ,(k)Δt，其中，Δfb,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力耦合的附加体积力的迭代误差；fb,(k)为第n+1个时间步下，第k个迭代步，与压力耦合的附加体积力；fb,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力耦合的附加体积力。
进一步地，所述步骤6)具体过程如下：
步骤6-1)将连续性条件表示为散度方程：
步骤6-2)将散度方程展开：则其中，Hn+1，(k+1)是离散动量方程中第n+1个时间步，第k+1个迭代步与压力pn+1相关的项。
进一步地，所述步骤7)中第n+1个时间步下、第k+1个迭代步的校正速度为un+1,(k+1)＝un+Gn+faΔt+fb,(k+1)Δt+Hn+1,(k+1)。
本发明所达到的有益效果：本方法解决了流场中存在浸入边界时流体流态计算的问题，改善了直接力法中存在的不足；充分考虑了附加体积力与压力的耦合，采用隐式迭代格式，计算时采用当前时间步中的压力，提高了浸入边界法数值模拟的准确性与可靠性；在附加体积力与压力迭代计算内进行速度校正，计算结果同时满足无滑边界条件与连续性条件，计算更符合实际。
附图说明
图1是圆柱绕流计算网格示意图；
图2是Re＝40时流场处于稳态时的流线图；
图3是Re＝100时流场处于稳态时的流线图；
图4是Re＝200时流场涡量图；
图5是Re＝1000时流场涡量图；
图6是Re＝100时瞬时压强等值线图；
图7是Re＝200时瞬时压强等值线图；
图8是Re＝100时圆柱迎流驻点附近速度矢量图；
图9是Re＝200时圆柱迎流驻点附近速度矢量图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
本发明提供了一种同时满足无滑边界条件和连续性条件的浸入边界隐式迭代求解方法，具体包括以下步骤：
步骤1)采用粘性不可压流N-S控制方程，将动量方程沿特征线展开进行时间离散，利用有限元法进行空间离散，在离散后的控制方程中加入附加体积力。具体步骤如下：
步骤1-1)采用粘性不可压流N-S方程，具体为：
动量方程：
连续性方程：
其中，是梯度算子，τ是粘性应力张量，μ是动力粘滞系数，u是流体速度，t是时间，ρ是流体密度，p是压力，fb是重力等体积力。
步骤1-2)：将动量方程沿特征线展开，得到时间离散后的动量方程，加入附加体积力，具体为：
其中，Δu为第n+1个和第n个时间步速度的差，Δt是计算时间步长，pn+1为第n+1个时间步的压力，fn+1为第n+1个时间步的附加体积力。
引入两种算子Gn和Hn+1，具体为：
其中，Gn为时间离散后的动量方程中，第n个时间步下，与速度un相关的项；Hn+1为时间离散后的动量方程中第n+1个时间步下，与压力pn+1相关的项。
则离散后的动量方程可表示为：
Δu＝Gn+Hn+1+fn+1Δt
即un+1＝un+Gn+Hn+1+fn+1Δt
其中，un+1为第n+1个时间步的速度。
步骤2)采用直接力法，使得流场插值速度等于边界上期望的速度分布，由此得到求解式。具体步骤如下：
步骤2-1)定义大写字母代表浸入边界点上的物理量，定义小写字母代表笛卡尔网格交界点上的物理量。
步骤2-2)定义两个函数，将定义在网格交界点上的物理量插值到浸入边界点上，D(Φ)将定义在浸入边界点上的物理量投射到网格交界点上，为定义在网格交界点上的物理量，Φ是定义在浸入边界点上的物理量。
其中，u(xi)为网格交界点xi上的速度，U(Xj)为浸入边界点Xj上的速度，δ(xi-Xj)为网格交界点xi与浸入边界点Xj上的物理量进行插值计算时所采用的插值函数，Ω为插值计算中涉及到的网格单元，ΔΩi是网格交界点xi的体积。
其中，F(Xj)为浸入边界点Xj上的力，f(xi)为xi网格交界点上的力，ΔVj是浸入边界点Xj的体积，N为浸入边界点个数。
步骤2-3)采用直接力法，使得流场插值速度等于边界上期望的速度分布，即满足无滑边界条件，浸入边界点上的速度等于笛卡尔网格交界点的插值速度：I(un+1)＝Ud，其中，Ud是第n+1个时间步边界上期望的速度分布。
则根据步骤1)中的动量离散方程，取算子：
I(un+Gn+Hn+1+fn+1Δt)＝Ud
即I(fn+1Δt)＝Ud-I(un+Gn+Hn+1)
进一步地，取算子：
D[I(fn+1Δt)]＝D[Ud-I(un+Gn+Hn+1)]
若D[I(fn+1Δt)]＝fn+1Δt，那么可以得到求解式：
fn+1Δt＝D[Ud-I(un+Gn+Hn+1)]
步骤3)根据步骤2)的求解式将附加体积力fn+1按是否与压力耦合分为两部分：fn+1Δt＝faΔt+fbΔt,其中，fa为第n+1个时间步中与压力不耦合的附加体积力，fb为第n+1个时间步中与压力耦合的附加体积力。
根据步骤2)中的求解式，则
faΔt＝D[Ud-I(un+Gn)]
fbΔt＝-D[I(Hn+1)]
步骤4)采用分步法求解格式，根据无滑边界条件，对与压力不耦合的附加体积力fa进行迭代计算，得到第n+1个时间步下的faΔt。具体步骤如下：
步骤4.1，假设fn,(k+1)Δt＝fn,(k)Δt+Δfn,(k)Δt
其中，fn,(k+1)为第n个时间步下，第k+1个迭代步的力，fn,(k)为第n个时间步下，第k个迭代步的力，Δfn,(k)为第n个时间步下，第k个和第k+1个迭代步力的迭代误差。
步骤4.2，根据无滑边界条件：
Δfa，(k)Δt＝D[Ud-I(un+Gn+fa，(k)Δt)]
fa,(k+1)Δt＝fa,(k)Δt+Δfa,(k)Δt
其中，Δfa,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力不耦合的附加体积力的迭代误差；fa,(k)为第n+1个时间步下，第k个迭代步，与压力不耦合的附加体积力；fa,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力不耦合的附加体积力。
迭代计算，直至收敛，得到第n+1个时间步faΔt的值。
步骤5)，根据无滑边界条件，对与压力耦合的附加体积力fb进行迭代求解。结合第k个迭代步的fb,(k)Δt和Hn+1,(k)，得到第k+1个迭代步的fb,(k+1)Δt。
Δfb，(k)Δt＝D[Ud-I(un+Gn+faΔt+fb，(k)Δt+Hn+1，(k))]
fb,(k+1)Δt＝fb,(k)Δt+Δfb,(k)Δt
其中，Δfb,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力耦合的附加体积力的迭代误差；fb,(k)为第n+1个时间步下，第k个迭代步，与压力耦合的附加体积力；fb,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力耦合的附加体积力。
步骤6)，结合步骤5)中第k+1个迭代步中的fb,(k+1)Δt，根据连续性条件，得到压力第k+1个迭代步的值Hn+1,(k+1)。具体步骤如下：
步骤6.1，连续性条件可以表示为散度方程：
步骤6.2，将散度方程展开：
则
其中，Hn+1,(k+1)是离散动量方程中第n+1个时间步，第k+1个迭代步与压力pn+1相关的项。
步骤7)根据步骤4)至步骤6)中的计算结果，得到第n+1个时间步下、第k+1个迭代步的校正速度un+1,(k+1)，公式为un+1,(k+1)＝un+Gn+faΔt+fb,(k+1)Δt+Hn+1,(k+1)；
步骤8)重复步骤4)至步骤7)，迭代计算直至收敛，得到第n+1个时间步的流速un+1；
步骤9)重复步骤4)至步骤8)，进行下一个时间步的迭代计算，直至计算完所有时间步。
实施例：以经典圆柱绕流算例，在流场中放置一个的圆柱，研究不同雷诺数下流体的流态。采用本发明提供的方法，为圆柱及流场建立单元模型，如图1所示，共划分流场单元215040个，浸入边界点418个。
图2和图3是雷诺数分别等于40和100，流场处于稳态时的流线图。图4和图5是雷诺数分别等于200和1000，产生卡门涡街时的涡量图。图6和图7是雷诺数分别等于100和200时，瞬时压强等值线图。图8和图9是雷诺数分别等于100和200时，圆柱迎流侧驻点附近的速度矢量图。基于本发明所提供的一种同时满足无滑边界条件和连续性条件的浸入边界隐式迭代求解方法模拟圆柱绕流过程，能清楚地描述固定圆柱处于不同雷诺数的流场中，流体流态的变化过程，可以很好的分析固体处于流场中对流场的影响，流场的速度、流线、涡量和压强分布等都能很直观的展示。
综上所述，采用基于浸入边界法的隐式迭代求解方法，能够实现存在浸入边界时流场流态的模拟。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。
Claims (9)
1.一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，包括如下步骤：
步骤1)采用粘性不可压流N-S控制方程，将动量方程沿特征线展开进行时间离散，利用有限元法进行空间离散，在离散后的控制方程中加入附加体积力；
步骤2)采用直接力法求解附加体积力，使得流场插值速度等于边界上期望的速度分布，得到求解式；
步骤3)根据步骤2)的求解式，将附加体积力按是否与压力耦合分为两部分，记为fa、fb，fa为第n+1个时间步中与压力不耦合的附加体积力，fb为第n+1个时间步中与压力耦合的附加体积力；
步骤4)采用分步法求解格式，根据无滑边界条件，对与压力不耦合的附加体积力fa进行迭代计算，得到第n+1个时间步下的fa；
步骤5)根据无滑边界条件，对与压力耦合的附加体积力fb进行迭代求解，得到第n个时间步下、第k+1个迭代步的fb,(k+1)；
步骤6)结合步骤5)中第k+1个迭代步中的fb,(k+1)，根据连续性条件，得到压力第n+1个时间步下、第k+1个迭代步的值Hn+1,(k+1)；
步骤7)根据步骤4)-6)中的计算结果，得到第n+1个时间步下、第k+1个迭代步的校正速度un+1,(k+1)；
步骤8)重复步骤4)-7)，迭代计算直至收敛，得到第n+1个时间步的流速un+1；
步骤9)重复步骤4)-8)，进行下一个时间步的迭代计算，直至计算完所有时间步。
2.根据权利要求1所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤1)具体步骤如下：
步骤1-1)采用粘性不可压流N-S控制方程，具体为：
动量方程：
连续性方程：▽·u＝0；
其中，▽为梯度算子，τ表示粘性应力张量，μ表示动力粘滞系数，u表示流体速度，t表示时间，ρ是流体密度，p是压力，fb为重力等体积力。
步骤1-2)将动量方程沿特征线展开，得到时间离散后的动量方程，加入附加体积力，具体为：
其中，Δu为第n+1个和第n个时间步速度的差，Δt是计算时间步长，pn+1为第n+1个时间步的压力，fn+1为第n+1个时间步的附加体积力。
3.根据权利要求2所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，引入两种算子Gn和Hn+1将离散后的动量方程重新表示，
其中，Gn为时间离散后的动量方程中，第n个时间步中，与速度un相关的项；Hn+1为时间离散后的动量方程中第n+1个时间步中，与压力pn+1相关的项；
离散后的动量方程可表示为：Δu＝Gn+Hn+1+fn+1Δt，即un+1＝un+Gn+Hn+1+fn+1Δt，其中，un +1为第n+1个时间步的速度。
4.根据权利要求3所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤2)具体步骤如下：
步骤2-1)定义大写字母代表浸入边界点上的物理量，定义小写字母代表笛卡尔网格交界点上的物理量；
步骤2-2)定义两个函数： 其中，I(φ)表示在网格交界点上的物理量插值到浸入边界点上，D(Φ)表示在浸入边界点上的物理量投射到网格交界点上，φ表示在网格交界点上的物理量，Φ表示在浸入边界点上的物理量，u(xi)为网格交界点xi上的速度，U(Xj)为浸入边界点Xj上的速度，δ(xi-Xj)为网格交界点xi与浸入边界点Xj上的物理量进行插值计算时所采用的插值函数，Ω为插值计算中涉及到的网格单元，ΔΩi是网格交界点xi的体积，F(Xj)为浸入边界点Xj上的力，f(xi)为网格交界点xi上的力，ΔVj是浸入边界点Xj的体积，N为浸入边界点个数；
步骤2-3)采用直接力法，使得流场插值速度等于边界上期望的速度分布，即满足无滑边界条件，浸入边界点上的速度等于笛卡尔网格交界点的插值速度：I(un+1)＝Ud，其中，Ud是第n+1个时间步边界上期望的速度分布；
根据动量离散方程，取算子：
I(un+Gn+Hn+1+fn+1Δt)＝Ud，即I(fn+1Δt)＝Ud-I(un+Gn+Hn+1)
D[I(fn+1Δt)]＝D[Ud-I(un+Gn+Hn+1)]
若D[I(fn+1Δt)]＝fn+1Δt，得到求解式fn+1Δt＝D[Ud-I(un+Gn+Hn+1)]。
5.根据权利要求4所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤3)中将附加体积力fn+1按是否与压力耦合分为两部分：fn+1Δt＝faΔt+fbΔt，其中，fa为第n+1个时间步中与压力不耦合的附加体积力，fb为第n+1个时间步中与压力耦合的附加体积力，faΔt＝D[Ud-I(un+Gn)]，fbΔt＝-D[I(Hn+1)]。
6.根据权利要求5所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤4)的具体步骤为：
步骤4-1)令fn,(k+1)Δt＝fn,(k)Δt+Δfn,(k)Δt，其中，fn,(k+1)为第n个时间步下，第k+1个迭代步的力，fn,(k)为第n个时间步下，第k个迭代步的力，Δfn,(k)为第n个时间步下，第k个和第k+1个迭代步力的迭代误差；
步骤4-2)根据无滑边界条件：Δfa,(k)Δt＝D[Ud-I(un+Gn+fa,(k)Δt)]，fa,(k+1)Δt＝fa,(k)Δt+Δfa,(k)Δt，其中，Δfa,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力不耦合的附加体积力的迭代误差；fa,(k)为第n+1个时间步下，第k个迭代步，与压力不耦合的附加体积力；fa,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力不耦合的附加体积力；
迭代计算，直至收敛，得到第n+1个时间步faΔt的值。
7.根据权利要求6所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤5)具体过程如下：
Δfb,(k)Δt＝D[Ud-I(un+Gn+faΔt+fb,(k)Δt+Hn+1,(k))]，fb,(k+1)Δt＝fb,(k)Δt+Δfb,(k)Δt，其中，Δfb,(k)为第n+1个时间步下，第k个迭代步和第k+1个迭代步中，与压力耦合的附加体积力的迭代误差；fb,(k)为第n+1个时间步下，第k个迭代步，与压力耦合的附加体积力；fb ,(k+1)为第n+1个时间步下，第k+1个迭代步，与压力耦合的附加体积力。
8.根据权利要求7所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤6)具体过程如下：
步骤6-1)将连续性条件表示为散度方程：▽·un+1＝0
步骤6-2)将散度方程展开：▽·[un+Gn+faΔt+fb,(k+1)Δt+Hn+1,(k+1)]＝0，则▽·(Hn +1,(k+1))＝-▽·[un+Gn+faΔt+fb,(k+1)Δt]，其中，Hn+1,(k+1)是离散动量方程中第n+1个时间步，第k+1个迭代步与压力pn+1相关的项。
9.根据权利要求8所述的一种满足无滑边界条件和连续性条件浸入边界隐式迭代求解法，其特征是，所述步骤7)中第n+1个时间步下、第k+1个迭代步的校正速度为un+1,(k+1)＝un+Gn+faΔt+fb,(k+1)Δt+Hn+1,(k+1)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710628403.4A CN107423511A (zh) | 2017-07-28 | 2017-07-28 | 满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710628403.4A CN107423511A (zh) | 2017-07-28 | 2017-07-28 | 满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107423511A true CN107423511A (zh) | 2017-12-01 |
Family
ID=60431606
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710628403.4A Pending CN107423511A (zh) | 2017-07-28 | 2017-07-28 | 满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107423511A (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108595759A (zh) * | 2018-03-22 | 2018-09-28 | 南京航空航天大学 | 一种基于OpenFOAM再开发的动边界问题计算方法 |
| CN109033712A (zh) * | 2018-08-31 | 2018-12-18 | 大连海事大学 | 一种基于物理的海上溢油实时交互模拟方法 |
| CN109446591A (zh) * | 2018-10-09 | 2019-03-08 | 苏州科技大学 | 考虑结构-载荷-边界耦合影响的结构参数优化设计方法 |
| CN111159853A (zh) * | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | 一种针对高雷诺数粘性流动问题的边界层高精度处理方法 |
| CN112131802A (zh) * | 2020-08-14 | 2020-12-25 | 山东大学 | 一种基于近场动力学的裂隙岩体渗流模拟方法及系统 |
| CN114117959A (zh) * | 2021-11-22 | 2022-03-01 | 北京航空航天大学 | 针对进发匹配的压气机稳定边界预测方法、设备及介质 |
| CN115577600A (zh) * | 2022-11-21 | 2023-01-06 | 广东电网有限责任公司中山供电局 | 一种单芯电终端头固-流两相对流换热计算方法和系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080063987A1 (en) * | 2006-09-13 | 2008-03-13 | Shinichi Ito | Immersion supporting plate cleaning method and a pattern forming method |
| CN101833605A (zh) * | 2010-04-29 | 2010-09-15 | 浙江工业大学 | 基于流体体积模型模具微细流道磨料流精密加工控制方法 |
| CN103970989A (zh) * | 2014-04-15 | 2014-08-06 | 昆明理工大学 | 一种基于流固界面一致条件的浸入边界流场计算方法 |
| CN105760602A (zh) * | 2015-12-30 | 2016-07-13 | 南京航空航天大学 | 一种有限体积加权基本无振荡格式的全流场数值模拟方法 |
-
2017
- 2017-07-28 CN CN201710628403.4A patent/CN107423511A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080063987A1 (en) * | 2006-09-13 | 2008-03-13 | Shinichi Ito | Immersion supporting plate cleaning method and a pattern forming method |
| CN101833605A (zh) * | 2010-04-29 | 2010-09-15 | 浙江工业大学 | 基于流体体积模型模具微细流道磨料流精密加工控制方法 |
| CN103970989A (zh) * | 2014-04-15 | 2014-08-06 | 昆明理工大学 | 一种基于流固界面一致条件的浸入边界流场计算方法 |
| CN105760602A (zh) * | 2015-12-30 | 2016-07-13 | 南京航空航天大学 | 一种有限体积加权基本无振荡格式的全流场数值模拟方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张修忠 等: "浅水方程数值计算方法的研究", 《水科学进展》 * |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108595759A (zh) * | 2018-03-22 | 2018-09-28 | 南京航空航天大学 | 一种基于OpenFOAM再开发的动边界问题计算方法 |
| CN109033712A (zh) * | 2018-08-31 | 2018-12-18 | 大连海事大学 | 一种基于物理的海上溢油实时交互模拟方法 |
| CN109033712B (zh) * | 2018-08-31 | 2020-05-29 | 大连海事大学 | 一种基于物理的海上溢油实时交互模拟方法 |
| CN109446591A (zh) * | 2018-10-09 | 2019-03-08 | 苏州科技大学 | 考虑结构-载荷-边界耦合影响的结构参数优化设计方法 |
| CN111159853A (zh) * | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | 一种针对高雷诺数粘性流动问题的边界层高精度处理方法 |
| CN112131802A (zh) * | 2020-08-14 | 2020-12-25 | 山东大学 | 一种基于近场动力学的裂隙岩体渗流模拟方法及系统 |
| CN112131802B (zh) * | 2020-08-14 | 2023-09-05 | 山东大学 | 一种基于近场动力学的裂隙岩体渗流模拟方法及系统 |
| CN114117959A (zh) * | 2021-11-22 | 2022-03-01 | 北京航空航天大学 | 针对进发匹配的压气机稳定边界预测方法、设备及介质 |
| CN115577600A (zh) * | 2022-11-21 | 2023-01-06 | 广东电网有限责任公司中山供电局 | 一种单芯电终端头固-流两相对流换热计算方法和系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107423511A (zh) | 满足无滑边界条件和连续性条件浸入边界隐式迭代求解法 | |
| Abadie et al. | On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks | |
| Akkerman et al. | Isogeometric analysis of free-surface flow | |
| CN107729691B (zh) | 一种稀薄连续统一的气体流动特性数值模拟方法 | |
| CN105975645B (zh) | 一种基于多步的含激波区域飞行器流场快速计算方法 | |
| Yang et al. | A strongly coupled, embedded-boundary method for fluid–structure interactions of elastically mounted rigid bodies | |
| Vierendeels et al. | Implicit coupling of partitioned fluid–structure interaction problems with reduced order models | |
| Francois et al. | Computations of drop dynamics with the immersed boundary method, part 1: numerical algorithm and buoyancy-induced effect | |
| Yuan et al. | An immersed-boundary method based on the gas kinetic BGK scheme for incompressible viscous flow | |
| Jaiman et al. | Transient fluid–structure interaction with non-matching spatial and temporal discretizations | |
| WO2012031398A1 (en) | Numerical method for simulating subsonic flows based on euler equations in lagrangian formulation | |
| Ii et al. | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach | |
| Jin et al. | A unified moving grid gas-kinetic method in Eulerian space for viscous flow computation | |
| Zhang et al. | A sharp interface method for SPH | |
| CN110532727A (zh) | 可用于常见非牛顿流体的数值模拟方法 | |
| Luo et al. | Curvature boundary condition for a moving contact line | |
| Zhang et al. | Development of a VOF+ LS+ SPP method based on FLUENT for simulating bubble behaviors in the electric field | |
| Brehm et al. | Towards a viscous wall model for immersed boundary methods | |
| Rozza et al. | Reduced basis methods and a posteriori error estimators for heat transfer problems | |
| JP7042562B2 (ja) | 非圧縮過渡および定常状態ナビエ－ストークス方程式のための最適圧力射影法 | |
| Ni et al. | Remapping-free ALE-type kinetic method for flow computations | |
| Aldlemy et al. | Adaptive mesh refinement immersed boundary method for simulations of laminar flows past a moving thin elastic structure | |
| Sai et al. | General purpose versus special algorithms for high‐speed flows with shocks | |
| JP3587827B2 (ja) | 翼形性能の推定方法および翼形性能の推定プログラム | |
| McIntyre et al. | The immersed boundary method for water entry simulation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20171201 |
|
| RJ01 | Rejection of invention patent application after publication |