改进的移动边界方案，用于晶格玻尔兹曼方法中新出现的流体节点

《Computers & Fluids》：Improved moving boundary schemes for emerging fluid nodes in the lattice Boltzmann method

　　平林信一郎
东京大学前沿科学学院海洋技术、政策与环境系，日本柏市柏之滨5-1-5

**摘要**
在格子玻尔兹曼方法中，估算从固态转变为液态的新生节点的分布函数仍然是一个关键难题。本文提出了一种改进的处理方法，通过对靠近固液界面的固态节点应用虚拟流体和虚拟边界条件来处理这些节点。为了验证所提方法的有效性，对受到强迫横向振荡影响的圆柱体和方形圆柱体周围的流体进行了二维模拟，并与传统方法进行了比较。结果表明，与传统方法相比，所提方法显著抑制了流体力中的虚假振荡。这些振荡源自新生流体节点处局部力的不连续性，而所提方法在缓解这些效应方面非常有效。在短期模拟中，所提方法估算的力与刚性网格参考结果几乎相同；然而，在长期模拟中，由于时间序列中的相位差异，差异可能会变得显著。敏感性分析表明，无论采用何种填充方案，这种相位差异都可能发生，这取决于流动特性的复杂性。

**1. 引言**
流体与运动固体之间的相互作用在广泛的比例尺度上都有应用，从悬浮液中细颗粒的微观运动到大型海上浮动结构的动力学。前者的典型雷诺数约为O(10–3)，后者则高达O(10^7)。虽然数值模拟是分析这些现象的有效工具，但仍存在重大的计算挑战。

首先，固液相互作用的特点是边界是移动的。一种简单的方法是使用与固体一起刚性移动的贴体网格，从而可以在相对参考框架中求解流体方程。然而，这种方法对于实际应用中常见的多体系统来说不切实际。在这种情况下，必须重新网格化以适应物体的运动。但这会显著增加计算成本，因为需要在大量网格节点上进行广泛插值。

另一种方法是使用固定网格并结合自适应边界条件来处理移动的固液界面。这种方法特别适用于具有简单解析描述的几何形状，如球体、圆柱体和方形杆，因为可以容易确定边界位置。由于在固定网格上计算的固有简单性，计算成本与固定物体情况相当，因为避免了复杂的全局插值。然而，边界条件的实现变得更加复杂，直接影响模拟的总体精度。因此，显然需要更稳健和准确的移动边界模型。

第二个挑战是将这些方法应用于高雷诺数问题，这些问题通常需要大量的计算资源。并行计算通过提高处理速度和可用内存提供了一个有前景的解决方案。然而，基于纳维-斯托克斯方程的传统计算流体动力学（CFD）需要求解整个计算域内的压力泊松方程。这一方程的全局性质往往阻碍了并行化的效率，随着系统规模的增加，并行化效率往往会降低。

格子玻尔兹曼方法（LBM）是一种能够克服这些挑战的替代数值框架。LBM通过粒子分布函数的演化来表示流体动力学。基于玻尔兹曼方程，分布函数流向相邻节点并在局部发生碰撞。由于在曲线或复杂界面上实现边界条件的固有简单性，LBM已成功应用于涉及复杂几何形状的各种问题，如多孔介质[1,2]和悬浮液[3,4]。此外，由于避免了压力泊松方程，该算法具有高度局部性，非常适合大规模并行计算[5,6]。

在LBM文献中，已经提出了几种处理移动边界的方法。Ladd[3]提出了一种在节点中间点将分布函数反弹的方案，应用于悬浮颗粒的运动。Filippova和H?nel[7]以及Bouzidi等人[8]随后通过采用基于插值的反弹方案来更准确地跟踪边界位置，进一步改进了这种方法。为了进一步提高精度，Kao和Yang[9]提出了一种无插值的处理方法。

尽管在确定固液边界附近的分布函数方面取得了这些进展，但关于处理相变节点的问题仍然存在关键挑战。当固体通过固定网格系统移动时，前沿附近的流体节点被固体吞噬，转变为固体节点；相反，尾缘的固体节点被固体离开，成为新的流体节点。这两种转变都会引入虚假波动，从而导致计算出的流体力出现误差。特别是，新生流体节点的分布函数是未知的，必须重新构建。这一过程常常引入显著的数值不准确性，尤其是在估算作用在固体上的流体力时。这种用于估算新生流体节点分布函数的“填充”算法已被广泛研究，并在Wen等[10]和Ginzburg[11]中引入。特别是Ginzburg[11]提出了一种用于可变形边界的填充算法。最简单的方法是使用来自周围流体节点的局部平均分布函数。或者，通常使用局部边界速度和局部平均密度构建平衡分布函数。Lallemand和Luo[12]提出的外推方法也是一种可行的替代方案。另一种方法是将分布函数分配给固体内部的节点，并通过与流体节点相同的碰撞和流动步骤进行更新[[13],[14],[15]]。Lallemand和Luo[12]报告说，这些不同的新生流体节点方案得到了可比的结果。然而，由于靠近边界的流体节点的分布函数直接决定了局部动量交换，填充方案的选择会显著影响计算出的流体力。当大量节点同时从固态转变为液态时，这种效应尤为明显，导致显著的虚假力波动。Tao等[16]比较了五种不同的“填充”算法，但其评估主要集中在力波动的大小上，而不是每种算法的总体精度。最近，He等[17]提出了一种结合平衡分布函数和非平衡分量加权的初始化方法。尽管他们的方法比传统方法具有优势，但其预测大位移运动时流体力的精度尚未得到充分验证。作为传统反弹方案的替代方案，Krithivasan等人[18]提出了一种扩散反弹方案，并证明使用这种新方案的填充算法可以减少力波动。然而，建议扩散反弹方案不能保证无滑移条件，并表现出某种程度的边界滑移[10]，这可能会改变流动特性。上述比较研究仅限于圆形或球形几何形状，没有涉及物体移动时新生节点数量显著增加的方形圆柱体。

移动边界的质量守恒也是一个关键问题，可能会影响力评估和流动重建。Yin等[19]指出，由于填充过程引起的局部质量泄漏是不可避免的，并建议不应将其包括在动量评估中，尽管其对全局流场的影响可能有限。最近，Xu等[20]从理论上证明，局部质量泄漏源于固液边界处壁切割链接的内在缺陷。这一见解对于移动边界问题尤为重要，因为边界位置相对于格子会随时间变化。

本研究对战绕受迫振荡的圆形和方形物体流动的LBM模拟中新生节点的分布函数重建方法进行了比较评估。特别是，通过将其与与物体同步移动的固定刚性网格配置进行基准测试，澄清了每种方法的精度和数值特性。此外，我还提出了一种新方法，并阐明了其各自的优点和局限性。

**2. 格子玻尔兹曼方法的简要回顾**
LBM的详细信息在文献中有大量记载[4,6]，这里仅简要回顾与本研究相关的基本方面。在本研究中，采用了使用Bhatnagar-Gross-Krook（BGK）近似的单松弛时间模型LBM。控制方程表示为：
(1)
fα(x+cαδt,t+δt) = fα(x,t) – 1/τρ(x,t) [fα(x,t) – fαeq(x,t)]
其中f是流体粒子的分布函数，feq是平衡分布函数，x是位置，t是时间，δt是时间步长，ρ是流体的质量密度，τ是松弛因子，α是流体粒子速度矢量的方向，cα是离散速度矢量。质量密度ρ和流体动量ρu（其中u是流体速度矢量）分别作为分布函数的零阶和一阶矩计算：
(2)
ρ = ∑αfα,
(3)
ρu = ∑αcαfα.
方程(1)通过以下两个步骤实现：
(4)
f?α(x,t) = fα(x,t) – 1/τ[fα(x,t) – fαeq(x,t)],
(5)
fα(x+cαδt,t+δt) = f?α(x,t).
方程(4)和(5)分别表示碰撞和流动步骤。f?α(x,t)表示碰撞后的分布函数，表示局部松弛后的流体状态。

**3. 固液边界处理的描述**
**3.1. 反弹方案**
在更新靠近固液界面的流体节点的分布函数时，从相邻固态节点流出的粒子数量是未知的。为了解决这个问题，Filippova和H?nel[7]（以下简称FH）提出了一种考虑在边界处反弹分布函数的方案。从固态节点流向流体节点的粒子数量通过以下线性插值确定：
(6)
fαˉ(xs,t) = (1–χ)f?α(xf,t) + χfα*(xs,t) + 2wαρ / (3c2cαˉ·ub),
其中上划线表示相反方向，wα是权重因子，cis是声速。xf表示靠近边界的流体节点的位置（以下简称边界流体节点），xs=xf+cαδt表示相邻的固态节点的位置（以下简称边界固态节点）。ub是固体物体的速度，f*是在固态节点定义的虚构平衡分布函数：
(7)
fα*(xs,t) = wαρ(xf,t) [1 + 3c2cα·usf + 9/2c?(cα·uf)2 – 3/2c2uf·uf].
其中uf是xf处的速度。usf和χ的组合定义如下：
(8)
usf = u(xf + λcαˉδt,t),
χ = (2d–1)τ^(-1) – λ (for d < 1/2),
(9)
usf = (d–1)duf + 1/d ub,
χ = (2d–1)τ (for d ≥ 1/2),
其中d表示从边界流体节点xf到固液边界的α方向上的分数距离，按格子间距归一化。FH在方程(8)中采用λ=0。Mei等[21]（以下简称MLS）指出，当τ趋近于1时，FH模型本质上是不稳定的。为了确保数值稳定性，他们提出λ=1。MLS还建议λ=2可以进一步提高稳定性，但实际上这是不必要的。

**3.2. 流体力的估算**
作用在固体上的流体力使用动量交换方法[22,23]进行评估。这种方法通过反弹过程中传递的动量来计算局部力：
(10)
F|xf| = ∑α[fα(xf,t) + fαˉ(xs,t)]cαδx / δt.
其中δx是格子间距。通过对所有边界流体节点xf的局部力求和，获得作用在固体上的总流体力。Kaushal等[24]提出了一种基于雷诺输运定理的更通用评估算法，并显示出比相同网格分辨率下的动量交换算法更好的精度。在本研究中，由于简单性，采用了动量交换算法，并确保了足够的网格分辨率以确保计算的有效性。

**4. 围绕振荡体的流动的数值模拟**
**4.1. 模型描述**
为了量化边界处理对流体力评估的影响，进行了围绕振荡体的LBM模拟。研究了两种不同的几何形状：圆柱体和方形圆柱体。当前的模拟采用了作者先前工作[25]中使用的相同计算参数和网格分辨率，当时已经建立了网格收敛性和数值有效性。水平域在x和y方向上分别为60D和40D，其中D是圆柱体直径和方形圆柱体边长的特征长度尺度。模拟采用三维十五速度（D3Q15）模型，并在z方向上使用单层来表示二维流场。为了保持二维性，在z方向上应用周期性边界条件。在入口边界施加x方向上的均匀流速U，其他水平边界设置为出流边界。通过U和D计算出的雷诺数（Re）固定为500。为了确保数值精度，在物体附近细化了网格。根据[25]中的网格收敛性研究，确定了D/δx=80的网格分辨率。有关模拟模型的详细描述和验证，请参考[25]。固体物体最初位于矩形计算域的中心。在流体-固体边界上，实施了MLS的反弹机制来满足边界条件。固体在横向流动方向上经历振荡加速，由以下方程定义：(11)d2/dt = -A*ω^2*cosω*t*，其中A表示振荡幅度，ω表示角频率，上标星号表示通过特征速度U和特征长度尺度D归一化的无量纲参数。无量纲参数设置为A*=0.1和ω*=π。使用反弹机制（方程（6）–（9））并取λ=1更新边界流体节点上的分布函数，符合MLS策略。

4.2. 补充方案
对于新出现的流体节点，使用表1中列出的七种补充方案来估计未知的分布函数，这些方案适用于圆形和方形圆柱体。这些方法的详细描述将在以下部分提供。此外，作为一个没有移动边界的基准测试，使用固定于物体的刚性网格进行了模拟。在这种设置中，振荡加速方向与方程（11）中的方向相反。这种刚性网格情况作为验证移动边界结果的参考。

表1. 振荡物体的计算案例
- 圆形圆柱体案例
- 方形圆柱体案例

补充方案：
C1S1：局部平均（LA）
C2S2：外推（EX）
C3S3：平衡分布函数（EQ）
C4S4：非平衡分布函数（NEQ）
C5S5：VF
C6S6：VF + VB
C7S7：混合（VF+VB & LA）

4.2.1. 局部平均（LA）
在这种方法中，通过计算相邻流体节点的函数的空间平均值来近似未知的分布函数，表示为：(12)fα LA(x,t) = 1/Nf * ∑βfα(x+cβδt,t)，其中Nf表示相邻流体节点的数量。当相邻节点是固体节点时，它从方程（12）的求和中排除。

4.2.2. 外推（EX）
遵循Lallemand和Luo [12]的方法，使用空间外推来重构未知的分布函数，表示为：(13)fα EX(x,t) = 3fα(x+cαδt,t) - 3fα(x+2cαδt,t) + fα(x+3cαδt,t)，当fα(x?cαδt,t)为固体节点时；否则，(14)fα EX(x,t) = f?α(x?cαδt,t?δt)。

4.2.3. 平衡分布函数（EQ）
如[12]中所指，具有局部边界速度的平衡分布函数被分配给未知的分布函数，表示为：(15)fα EQ(x,t) = fα eq(ub,ρ0)，其中ρ0是水的标准密度。

4.2.4. 非平衡分布函数（NEQ）
遵循He等人[17]的方法，未知分布被重构为平衡分量和局部平均非平衡分量的总和：(16)fα NEQ(x,t) = fα eq(ub,ρ0) + 1/Nf * ∑βfαneq(x+cβδt,t)，其中fαneq定义为：(17)fαneq(x,t) = fα(x,t) - fα eq(x,t)。

4.2.5. 虚拟流体（VF）方案
在这种方法中，假设固体内部充满了“虚拟流体”，并且也将分布函数分配给固体节点[[13], [14], [15]]。此后，这种方法被称为虚拟流体（VF）方案。VF节点的分布函数在碰撞和流动步骤中的更新与流体节点相同。在流动步骤中，任何从固体节点流到流体节点的分布函数都被替换为第3.1节中显示的反弹机制得到的函数。图1(a)示意性地展示了流体节点和VF节点之间的分布函数流动。相反，从流体节点流到固体节点的分布函数用于更新内部固体节点的函数。因此，这些VF节点内的计算主要用于重构新出现节点的分布函数，并不直接影响固体外部流场。

下载：下载高分辨率图像（390KB）
下载：下载全尺寸图像

图1. (a) 虚拟流体（VF）方案和(b) 虚拟流体和虚拟边界（VF+VB）方案中的流动和反弹（BB）过程示意图。粗实线表示流体-固体边界，虚线表示虚拟边界（VBs）。方形表示流体节点；圆形表示固体节点。特别注意，空心方形和空心圆形分别表示边界流体节点（xf）和边界固体节点（xs）。请注意，对角线流动分量被省略，图中仅显示选定的分量。

当VF节点转变为流体节点时，它已经充满了分布函数。然而，Ginzburg [14]指出，使用从VF节点向固体边界的填充会导致数值波动，并建议在填充过程中用其他方法替换这些波动。这个问题将在4.2.6节中讨论。

4.2.6. 组合虚拟流体和虚拟边界（VF+VB）方案（新方法）
正如Ginzburg [14]针对VF方案指出的，从VF节点向流体-固体边界的填充是不现实的。这可以解释如下：当VF节点刚刚转变为流体节点时，它本质上位于流体-固体边界附近。虽然这些节点的流体力通常通过动量交换算法（方程（10）来估计，但VF方案中的分布函数没有考虑局部边界约束。这种缺乏约束可能导致朝向流体-固体边界的速度分量被高估。尽管流体-固体边界的效应会通过反弹过程随时间反映出来，但在转换后立即产生的不切实际的高估速度会导致冲击压力峰值。为了将这些边界效应纳入VF方案，在固体内部引入了一个虚拟边界（VB），如图1(b)所示。这个VB位于边界固体节点的物理边界对面的一侧，即内部。换句话说，边界固体节点在几何上位于物理边界和虚拟边界之间。在这些节点上，仅对VB应用反弹条件，而不対物理边界应用。通过在VB处应用反弹条件，可以在边界固体节点转变为流体节点之前有效地捕获局部边界约束。采用这种处理方式，边界固体节点在物理上能够良好地转变为边界流体节点。

这种方案（此后称为VF+VB方案）应用于边界固体节点。虽然内部的固体节点可以被视为普通的VF节点，但它们的特定处理并不重要，因为它们不立即经历相变，也不会影响边界固体节点。边界固体节点处的分布函数的数学公式与用于边界流体节点的公式相同。它包括碰撞和流动步骤（方程（4）和（5），并结合了反弹机制（方程（6）–（9）），将边界固体节点视为“虚拟”边界流体节点。由于VB的确切位置不是唯一确定的，在模拟中从固体节点到VB的归一化分数距离dVB可以取任意值。然而，考虑到VB在固体节点转变为流体节点后立即转变为物理边界——此时归一化分数距离d接近零——将dVB设置为零以保持这种转换过程中的一致性是合理的。因此，在本研究中，dVB被设置为零。严格来说，刚刚转变为流体节点的d在数值模拟中并不完全为零，这种差异可能会导致力评估中的微小数值峰值。然而，只要边界在一个时间步长内的移动与晶格间距相比足够小，这些效应就可以忽略不计。

这种VF+VB方案是对现有VF方案的改进，它调整了边界固体节点的填充情况，以更符合新的流体节点的实际条件。这种新方法在本研究中提出，并用于表1中的C6和S6案例。

4.2.7. VF+VB方案与局部平均（混合）结合
虽然VF+VB方案预期可以纠正朝向边界流体节点的分布函数，但剩余部分仍然依赖于VF方案。由于VF节点仅从流体节点接收单向流动，它们的分布函数的方向组成可能与周围流场不一致。这可能导致刚刚转变的流体节点上的速度不切实际。为了补偿这个问题，提出了一种混合方法。该方法的思路是利用VF+VB方案对动量大小的校正，同时保持速度与周围流场的方向一致性。在这种方法中，VF+VB方案用于密度（和压力）计算，同时对分布函数应用局部平均技术，表示为：(18)ρVF+VB = ∑αfαVF+VB，(19)ρLA = ∑αfαLA，(20)fαhybrid(xs,t) = ρVF+VB * ρLA * fαLA(xs,t)。对于局部平均，采用了与4.1.1相同的方案。这种混合方法应用于C7和S7案例。

5. 结果和讨论
5.1. 初始时期（0≤Ut/D≤4）
在本节中，讨论了初始时期0≤Ut/D≤4的结果。

图2显示了振荡圆形圆柱体的阻力系数的时间序列，并与固定于物体的刚性网格基准结果（以下简称基准）进行了比较。虽然所有案例都遵循与基准相同的总体趋势，但波动的幅度有所不同。C1（LA）、C2（EX）和C3（EQ）的波动最为明显，而C4（NEQ）、C5（VF）和C7（VF+VB & LA）的波动相对较小。相比之下，C6（VF+VB）的波动几乎可以忽略不计。同样，图3展示了升力系数的时间序列。C1的波动最为显著，其次是C2、C3、C5和C7，它们的波动水平相当。与阻力系数一样，C6中的波动明显得到了抑制。比较C3和C4案例，NEQ方案明显比EQ方案更能抑制力的波动。

下载：下载高分辨率图像（403KB）
下载：下载全尺寸图像

图2. 振荡圆形圆柱体的阻力系数的时间序列。虚线表示案例（a）C1、（b）C2、（c）C3、（d）C4、（e）C5、（f）C6和（g）C7，与刚性网格结果（实线）进行比较。

图3. 振荡圆形圆柱体的升力系数的时间序列。虚线表示案例（a）C1、（b）C2、（c）C3、（d）C4、（e）C5、（f）C6和（g）C7，与刚性网格结果（实线）进行比较。

图4显示了振荡方形圆柱体的阻力系数的时间序列，并与基准结果进行了比较。由于移动边界处理的不连续性，图中出现了较大的虚假峰值。与圆形物体相比，方形圆柱体的峰值幅度更大。这是因为在当前模拟中，方形圆柱体的边界始终与网格排列平行，而且由于方形物体的移动，单个时间步骤内新转换的流体节点数量更多。在S1和S2案例中观察到了显著的虚假峰值，而在S3、S5和S7案例中，这些波动相对较小。值得注意的是，在S4和S6案例中，峰值大幅减少。图5显示了升力系数的时间序列。与阻力系数类似，S1和S2案例中的峰值幅度较大，S3和S4案例中幅度适中，而在S5、S6和S7案例中进一步减少。图5显示，S6案例中的波动减少最为显著。

下载：下载高分辨率图像（633KB）
下载：下载全尺寸图像

图4. 振荡方形圆柱体的阻力系数的时间序列。虚线表示案例（a）S1、（b）S2、（c）S3、（d）S4、（e）S5、（f）S6和（g）S7，与刚性网格结果（实线）进行比较。

图5. 振荡方形圆柱体的升力系数的时间序列。虚线表示案例（a）S1、（b）S2、（c）S3、（d）S4、（e）S5、（f）S6和（g）S7，与刚性网格结果（实线）进行比较。

5.2.延长期（100≤Ut/D≤120）图6和图7展示了在100≤Ut/D≤120延长期间圆柱体的阻力和升力系数。虽然每种情况下的波动幅度与图2和图3中呈现的结果一致，但在案例C2、C3、C4、C5和C6中与基准值存在显著差异。在这些案例中，C2、C4、C5和C6的趋势与基准值相似，尽管它们的相位有所偏移。由于这些案例在图2和图3的初始阶段显示出与基准值相同的相位，随后的偏差表明当前方法在捕捉那些最初可以忽略但随时间累积变得显著的物理效应方面存在局限性。这可能是由于在重新填充过程中对分布函数的方向特性表示不准确造成的。例如，尽管由非平衡分布函数表示的粘性效应可能不是瞬时力的主要贡献因素，但它们通过改变流动模式（如涡脱落）来影响时间序列的相位。下载：下载高分辨率图像（473KB）下载：下载全尺寸图像图6. 在100≤Ut/D≤120区间内，振荡圆柱体的阻力系数时间序列。虚线代表案例(a) C1、(b) C2、(c) C3、(d) C4、(e) C5、(f) C6、(g) C7，与刚性网格结果（实线）进行比较。下载：下载高分辨率图像（606KB）下载：下载全尺寸图像图7. 在100≤Ut/D≤120区间内，振荡圆柱体的升力系数时间序列。虚线代表案例(a) C1、(b) C2、(c) C3、(d) C4、(e) C5、(f) C6、(g) C7，与刚性网格结果（实线）进行比较。相比之下，案例C1和C7准确跟随了参考趋势。这表明分布函数的平均处理在捕捉固有流动特性方面起着关键作用，确保了它们在长期积分过程中的保持。特别是，案例C7成功地重建了流体力，且波动最小。图8和图9显示了在Ut/D=100时的密度分布——即标准化压力。与图8中的基准结果相比，很明显案例C1和C7最准确地捕捉了流动趋势。案例C3中的流场与参考值相当接近，这与图6和图7中观察到的力量评估结果的合理一致性一致。相比之下，其他案例对涡脱落相关的压力振荡要么高估要么低估，这可能导致力量评估中的相位不准确。然而，如第5.4节后面讨论的，这些长期模拟中的趋势并不一定在所有条件下都成立，因为它们高度依赖于特定的流动模式。下载：下载高分辨率图像（426KB）下载：下载全尺寸图像图8. 在Ut/D=100时刚性网格情况下的密度（标准化压力）等高线图。下载：下载高分辨率图像（2MB）下载：下载全尺寸图像图9. 在Ut/D=100时，案例(a) C1、(b) C2、(c) C3、(d) C4、(e) C5、(f) C6和(g) C7的密度（标准化压力）等高线图。图10和图11展示了在100≤Ut/D≤120延长期间方形圆柱体的阻力和升力系数。与圆柱体案例不同，所有方法都准确跟随了总体趋势。这可能是因为方形圆柱体的涡脱落机制比圆柱体简单，因为其分离点固定，从而能够进行准确的力量评估。值得注意的是，在案例S3、S4、S5、S6和S7中，波动幅度显著减小，其中S6的波动最小。下载：下载高分辨率图像（826KB）下载：下载全尺寸图像图10. 在100≤Ut/D≤120区间内，振荡方形圆柱体的阻力系数时间序列。虚线代表案例(a) S1、(b) S2、(c) S3、(d) S4、(e) S5、(f) S6和(g) S7，与刚性网格结果（实线）进行比较。下载：下载高分辨率图像（806KB）下载：下载全尺寸图像图11. 在100≤Ut/D≤120区间内，振荡方形圆柱体的升力系数时间序列。虚线代表案例(a) S1、(b) S2、(c) S3、(d) S4、(e) S5、(f) S6和(g) S7，与刚性网格结果（实线）进行比较。5.3. 误差分析图12和图13分别展示了相对于基准值的圆柱体和方形圆柱体的均方根误差（RMSE）。请注意，在此上下文中，RMSE是通过对每个模拟的瞬时误差在特定时间窗口内积分到刚性网格基准值得到的单个标量。如图时间序列所示，由于VF+VB方案的有效波动抑制作用，在初始阶段（黑色条形）案例C6和S6显示出最低误差。然而，随着长期模拟的进行，案例C6的误差（白色条形）逐渐增加。另一方面，由于流动模式较为简单，案例S6在延长期间仍保持了其稳健的性能，如前面的图11所讨论的。下载：下载高分辨率图像（146KB）下载：下载全尺寸图像图12. 圆柱体所有案例的均方根误差（RMSE）：(a) 阻力系数和(b) 升力系数。黑色条形代表初始阶段（0≤Ut/D≤4），白色条形表示延长阶段（100≤Ut/D≤120）。下载：下载高分辨率图像（149KB）下载：下载全尺寸图像图13. 方形圆柱体所有案例的均方根误差（RMSE）：(a) 阻力系数和(b) 升力系数。黑色条形代表初始阶段（0≤Ut/D≤4），白色条形表示延长阶段（100≤Ut/D≤120）。5.4. 敏感性分析在本节中进行了敏感性分析，以验证当前研究结果的普遍性和稳健性。为了简化，目标仅限于圆形配置，并且只研究了案例C1、C4、C6和C7。5.4.1. 网格模型尽管本研究采用了D3Q15模型来离散化速度，但对其依赖于网格模型的问题存在担忧。实际上，Kandai等人（1999年）[26]指出，D3Q15可能会在3D模拟中产生非物理性的棋盘格波动，这可以通过使用D3Q19来抑制。其他研究表明，D3Q15和D3Q19模型的旋转不变性分别为差和中等，而D3Q27模型满足这一原则[27,28]。图14展示了D3Q15和D3Q19模型的阻力和升力系数相对于刚性网格基准值的RMSE比较。发现这两种模型的结果几乎相同，差异可以忽略不计。由于当前问题是二维的，先前研究中提出的问题并不明显。因此，我们可以得出结论，对于当前问题选择D3Q15模型是合理的。下载：下载高分辨率图像（168KB）下载：下载全尺寸图像图14. 案例C1、C4、C6和C7的圆柱体RMSE：(a) 阻力系数和(b) 升力系数。D3Q15模型的结果显示为黑色（初始阶段，0≤Ut/D≤4），D3Q19模型的结果显示为灰色和对角线阴影条形。5.4.2. 到VB的距离dVB在当前研究中，将到VB的标准化距离dVB设置为零，以保持VB过渡到物理边界时的一致性，如第4.2.6节所解释的。图15显示了案例C6（VF+VB方案）的阻力和升力系数的RMSE比较。在初始阶段，dVB=0的情况在阻力和升力系数方面都表现出最佳性能，这是预期的。dVB=0.5和1的结果几乎具有相同的RMSE，表明在这个范围内对dVB的敏感性已经饱和。下载：下载高分辨率图像（141KB）下载：下载全尺寸图像图15. 案例C6的圆柱体RMSE：(a) 阻力系数和(b) 升力系数。黑色条形代表初始阶段（0≤Ut/D≤4），白色条形表示延长阶段（100≤Ut/D≤120）。然而，在延长阶段，所有RMSE几乎相同，无论dVB的值如何，没有明显的依赖性。这些较大的RMSE主要是由时间序列相对于基准值的相位差异造成的。因此，dVB的影响相对可以忽略。图16展示了初期（0≤Ut/D≤4）情况下，具有不同dVB值的案例C6（VF+VB方案）和案例C5（VF方案）的阻力系数时间序列。可以看出，较小的dVB产生了较少的虚假峰值，这与图15中的结果一致。显然，VF+VB方案将边界固体节点上的分布函数校正到更符合物理现实的状态，从而抑制了重新填充过程中的非物理压力（或密度）跳动。因此，所提出的VF+VB方案预计可以缓解传统填充算法[19]中固有的局部质量泄漏问题。下载：下载高分辨率图像（222KB）下载：下载全尺寸图像图16. 振荡圆柱体C6（VF+VB方案）和C5（VF方案）的阻力系数时间序列。虚线代表案例(a) C6 with dVB=0、(b) C6 with dVB=0.25、(c) C6 with dVB=0.5 和 (d) C5，与刚性网格结果（实线）进行比较。5.4.3. 雷诺数图17展示了在Re = 100和500条件下，相对于刚性网格基准值的圆柱体和方形圆柱体的阻力和升力系数的RMSE比较。在Re = 100条件下，所有案例的RMSE都小于Re = 500条件下的RMSE。尽管RMSE的幅度有所不同，但在初始阶段，每种边界方案的相对特性是相同的，不受雷诺数的影响。下载：下载高分辨率图像（152KB）下载：下载全尺寸图像图17. 案例C1、C4、C6和C7的圆柱体RMSE：(a) 阻力系数和(b) 升力系数。Re = 500的结果显示为黑色（初始阶段，0≤Ut/D≤4），Re = 100的结果显示为灰色和对角线阴影条形。此外，即使在Re = 100条件下，案例C4和C6也在延长期间表现出良好的性能。可以推断，对于没有小涡流的相对简单流动，移动边界方法可以成功再现流场，并且VF+VB方案在长时间内保持了最佳的准确性。5.4.4. 振动参数通过比较三个振幅（A*=0.05,0.1,0.2）和四个频率（ω*/2π=0.1,0.25,0.5,1）的案例，研究了振动参数的影响。图18展示了相对于刚性网格基准值，在初始阶段（0≤Ut/D≤4）的阻力和升力系数的RMSE比较。随着振动幅度的增大，阻力和升力的RMSE分别显示出增加和减少的趋势。这可能与力的大小有关。然而，如图19所示，在延长阶段（100≤Ut/D≤120），这种一致性丧失了。较大的RMSE主要是由于每种方案与基准值之间的时间序列相位差距造成的。图20展示了在初始阶段（0≤Ut/D≤4）对于频率比较的阻力和升力系数的RMSE。与图18类似，可以看到每种边界方案的相对特性的一致性，表明了这一时期的研究结果的稳健性。然而，在延长阶段（100≤Ut/D≤120），如图21所示，一致性丧失，每种方案在不同频率条件下都可能显示出较大的RMSE。值得注意的是，在ω*/2π=0.1,0.25时，所有案例C1、C4、C6和C7的RMSE都较大。考虑到这些频率接近涡脱落频率，并且在接近共振条件下流动模式变得复杂，因此再现实时间序列的准确性可能取决于流动模式的复杂性。下载：下载高分辨率图像（149KB）下载：下载全尺寸图像图20. 在初始时间段（0≤Ut/D≤4）内，圆柱体C1、C4、C6和C7情况下的均方根误差（RMSE）：(a) 阻力系数和 (b) 升力系数。灰色、对角虚线、黑色和白色条分别对应 ω*/(2π)=0.1、0.25、0.5、1。下载：下载高分辨率图像（192KB）下载：下载全尺寸图像图21. 在延长时间段（100≤Ut/D≤120）内，圆柱体C1、C4、C6和C7情况下的均方根误差（RMSE）：(a) 阻力系数和 (b) 升力系数。灰色、对角虚线、黑色和白色条分别对应 ω*/(2π)=0.1、0.25、0.5、1。5.5. VF+VB方案的总体趋势在所有模拟中，VF+VB方案有效地抑制了力评估中的虚假峰值。特别是在短期模拟中，VF+VB方案对于圆形和方形配置都实现了足够低的RMSE。这一特性在各种流动条件下都表现得非常稳健，如第5.4节所示。然而，在长期模拟中，圆柱体情况由于明显的相位偏移而偏离了基准值，而方形圆柱体情况则与基准值保持一致。这两种配置的比较表明，重建精度高度依赖于流动特性。特别是，圆柱体的涡脱落对局部压力和剪切应力非常敏感。这使得在长期模拟中难以保持相位一致性，因为分离点并非像方形圆柱体那样是几何固定的。这一点也得到了第5.4节敏感性分析的支持，在该分析中，低雷诺数情况下的长期模拟保持了较低的RMSE值，因为此时涡脱落的影响相对较小。VF+VB与局部平均的混合方案在测试的填充算法中通常表现出中等的性能，这主要是由于VF+VB方案的特性。尽管在特定条件下，它在长期模拟中的性能优于VF+VB方案，但其性能仍然非常依赖于流动特性，这一点通过敏感性分析得到了证实。本研究提出的最优方案的定性总结见表2。表2. 本研究提出的最优填充方案。流动特性短期长期简单（方形，低Re等）VF+VBVF+VB复杂（圆形）VF+VBNA（特定情况）6. 结论本研究评估了现有的填充方法，并提出了一种处理LBM中移动流体-固体边界的新方法。使用刚性网格基准，评估了每种方法的数值准确性和虚假波动的幅度。结果表明，与传统的方法相比，所提出的VF+VB方案显著抑制了虚假波动，从而实现了高效的流动重建和力评估，适用于短期模拟。然而，在复杂流动模式下，由于局部小误差的累积，再现实象在延长时间段内往往会偏离基准值。所提出的方案为短期流动重建以及没有复杂涡脱落的情况（如方形圆柱体配置或低Re流动）提供了一个稳健的解决方案。由于本研究仅限于二维问题，因此将其扩展到三维情况是下一个挑战。CRediT作者贡献声明平林信一郎：写作 – 审稿与编辑，写作 – 原始草稿，可视化，验证，软件，方法论，研究，资金获取，正式分析，概念化。