一种用于加速边界条件约束浸没边界法（IBEM）的算法策略，适用于挤压壁几何体

《Computers & Fluids》：An algorithmic strategy to accelerate boundary-condition-enforced immersed boundary method for extruded wall geometries

　　摘要
在本研究中，提出了一种新颖的简化公式，用于计算边界条件强制扩散界面浸入边界方法（IBM）中校正算子的值，该方法适用于挤压壁几何形状。该简化公式利用了挤压壁几何形状在横向方向上的对称性，得出的解与原始的边界条件强制扩散界面IBM算法相同。通过对线性方程进行三步简化，将矩阵的大小从N×N降低到(N/r)×(N/r)，其中N表示拉格朗日点的数量，r表示横向方向的行数。通过多次计算证明，所提出的算法在计算效率上显著优于原始算法：校正算子的计算效率提高了2800倍，而整个IBM算法的计算效率提高了160倍，同时保持了计算精度。特别是，横向拉格朗日点数量的增加并不会导致计算成本的指数级上升，这是与其他高精度算法相比的一个独特特点。这种效率的提高显著扩展了边界条件强制扩散界面IBM在需要大量拉格朗日点和完全三维配置的高雷诺数流动模拟中的应用范围。

1. 引言
扩散界面浸入边界方法（IBM）[1],[2],[3]被广泛应用于计算流体动力学（CFD）中，用于研究围绕几何形状复杂且移动物体的流场。由于其相对稳定，并且可以通过考虑由一组拉格朗日点表示的壁面来轻松实现，因此得到了广泛应用，而无需生成复杂的计算网格。然而，扩散界面IBM的一个固有瓶颈是边界层的分辨率较低，因为外力通常通过使用狄拉克δ函数[4]在壁面周围扩散。为了克服这一问题，提出了几种IBM算法来在壁面处强制执行无滑移边界条件。Wu等人[5]提出了一种通过在壁面处求解代数方程来实现速度校正的方法。多方向强迫方法[6],[7]也提出了一种方法，在一个时间步长内施加多个外力，以逐步获得壁面处的期望速度。
在这方面，Pinelli等人[8]提出了边界条件强制IBM，该方法通过施加一个扩展算子校正来隐式确保扩散界面IBM外力的无滑移边界条件。与其他方法相比，它的优势在于壁面周围的速度场对算法收敛性的影响较小。尽管该校正成功地实现了壁面拉格朗日点处的边界条件，但求解每个拉格朗日点所需的校正算子的计算成本较高。随着拉格朗日点数量N的增加，计算成本也会增加。在需要大量拉格朗日点的大规模模拟中[9],[10]，校正算子的计算成本会迅速增加到流场计算的数十倍。先前的研究已经成功利用GPU加速了计算[11],[12]。然而，这种方法无法在纯CPU-based计算机上实现。此外，这种方法对于模拟高雷诺数流动时需要更多拉格朗日点的情况缺乏可扩展性，因为它没有解决可扩展性问题。
本文提出了一种利用简化计算来求解校正算子的算法方法，当壁面具有挤压几何形状时，该方法可以显著降低计算成本，同时保持计算精度。其主要创新之处在于通过算法简化来求解边界条件强制IBM中的校正算子，而不是依赖于硬件加速。因此，所提出的算法既可以在基于CPU的计算机上实现，也可以在基于GPU的计算机上实现，并且适用于表示壁面的更多拉格朗日点。此外，与现有的加速策略（如共轭梯度方法）相比，当壁面具有挤压几何形状时，所提出的方法具有更高的效率。该算法可以应用于满足以下条件的流体流动数值模拟：(a) 壁面具有挤压几何形状且保持相同的横截面几何形状；(b) 使用欧拉网格，该网格要么是均匀的，要么其沿挤压方向的拉伸率是恒定的。该算法也适用于部分具有(a)和(b)条件的壁面。在本文中，讨论了将边界条件强制IBM应用于不可压缩流动的情况。

2. 数值方法
2.1. 控制方程
不可压缩流动的质量守恒和动量守恒的控制方程写为：
(1) ??u = 0
(2) ρ?u?t + ρ??uu = ??P + ??τ + f
其中ρ是密度，u是速度，P是压力，τ是粘性应力张量，f是使用扩散界面浸入边界方法施加的假想外力[1]。非弹性壁面由一组拉格朗日点表示，通过在这些点上插值，可以通过施加假想外力f来修改欧拉网格上的物理量，从而使速度分布满足壁面的无滑移边界条件。本研究采用的数值方案基于一个内部热流分析代码FK3（例如[13],[14]），适用于不可压缩流动。

2.2. 传统的边界条件强制IBM
2.2.1. 基于力校正的边界条件强制IBM
在使用扩散界面IBM时，壁面由一组拉格朗日点表示，通过这些点，通过对欧拉网格上的物理量施加假想外力来修改它们，从而使所有物理量（如速度和压力）满足壁面的边界条件。假想外力的推导过程如下：
首先，在不考虑壁面存在的情况下求解动量守恒方程。然后，从欧拉网格的时间解中插值出壁面上拉格朗日点的速度：
(3) Us = ∑j∈Dsu?jδh|xj?Xs|Δxjδh|yj?Ys|Δyjδh|zj?Zs|Δzj
其中u?是考虑壁面存在之前的速度时间解，(Xs, Ys, Zs)分别是第s个拉格朗日点的速度和坐标，(xj,yj,zj), Δx, Δy, Δz分别是欧拉网格的坐标和网格大小。此外，Ds是欧拉网格中第s个拉格朗日点周围的小范围点集，狄拉克δ函数δh(r)表示为：
(4) δh(r) = 1
1 ? 3r2 + 1 (0 ≤ r ≤ 0.5)
16(5 ? 3r ? 1) (0.5 < r ≤ 1.5)
0otherwise
利用这个插值速度和每个拉格朗日点处移动壁面的速度Uw，可以使用以下公式计算拉格朗日点处的假想外力Fs：
(5) Fs = ρUw ? UsΔt
计算出Fs后，将其插值回欧拉网格：
(6) fs = ∑s∈DjFsδh|xj?Xs|Δxjδh|yj?Ys|Δyjδh|zj?Zs|ΔzjΔss?s
其中Δss和?s是连接Xs?1/2和Xs+1/2的弧长，校正算子由欧拉网格和拉格朗日点之间的位置关系确定，具体如下：
(7) Akl = ∑j∈Dlδh|xj?Xk|Δxjδh|xj?Xl|Δxj
(8) A? = 1
这里，矩阵A是一个大小为N×N的方阵，N是拉格朗日点的数量，Akl是矩阵A中第k行第l列的元素，表示拉格朗日点k和l之间的力干扰。这种校正隐式确保了壁面的无滑移边界条件，最初由Pinelli等人[8]提出。这种校正可以防止流体在壁面的意外渗透。最后，将插值后的外力fs施加到欧拉网格上：
(9) u = u? + fΔtρ
这种基于力校正的边界条件强制IBM的优势在于?的收敛性与周围的流场无关，因此整个模拟期间的计算成本是恒定的。
2.2.2. 基于速度校正的边界条件强制IBM
上述基于边界的IBM本质等同于Wu等人[5]提出的隐式基于速度校正的IBM，它校正每个拉格朗日点周围的速度分布，而不是校正外力。边界条件强制IBM通过速度校正的目的是使拉格朗日点处的边界速度与欧拉网格中的插值流体速度相匹配。目标速度分布可以表示为：
(10) Uw(Xk) = ∑j∈Dsujδh|xj?Xk|Δxj
其中Uw和u分别是拉格朗日点处的速度和未知的速度分布。未知的u可以表示为：
(11) u = u? + δu
其中δu是速度校正，对应于方程(9)中的外力。将u代入方程(10)得到：
(13) Uw(Xk) = ∑j∈Dku?jδh|xj?Xk|Δxj + ∑j∈Dkδujδh|xj?Xk|Δxj
此外，将方程(5)、(12)代入方程(6)并消除不相关的?s（仅在基于力校正的边界条件强制IBM中需要）得到：
(14) δuj = ∑s∈DsδUsδh|xj?Xs|ΔxjΔss
方程(13)、(14)现在得到：
(15) Uw(Xk) = ∑j∈Dsu?jδh|xj?Xk|Δxj + ∑s∈Dk∑j∈DsδUsδh|xj?Xk|Δxjδh|xj?Xk|ΔxjΔss
这里，δUs是未知变量，因此方程(15)可以重新定义为以下系统：
(16) AY = b
其中
(17) Y = δU1Δs1δU2Δs2?δUNΔsN
(18) b = Uw1Uw2?UwN?Du?1u?2?u?M
(19) D = δh|x1?X1|Δx1δh|x2?X1|Δx2…δh|xM?X1|ΔxMδh|x1?X2|Δx1?????δh|x1?XN|Δx1?…δh|xM?XN|ΔxM
方程(16)中的A相当于方程(7)中的A。然后可以利用获得的Y以以下形式将外力施加到流体上：
(20) f = ρδuΔt

2.3. 基于共轭梯度技术的IBM
前面提到的两种边界条件强制IBM都存在高计算成本的问题，因为在每个时间步长中，每次迭代都需要对A∈RN×N进行繁琐的矩阵评估。随着拉格朗日点数量(N)的增加，这种评估的计算成本会迅速增加到O(N2)，这使得边界条件强制IBM难以应用于高雷诺数流动或需要高分辨率的复杂几何形状周围的流动。可以通过避免在每次迭代中显式评估方程(8)或(16)来减轻计算负担。在Zhao等人[15]提出的基于共轭梯度技术的IBM中，通过仅应用密向量内积和稀疏矩阵-密向量乘法（SpMV），避免了显式评估方程(16)，每次迭代只需要LN次浮点运算（FLOPs）。这里，N和L分别表示拉格朗日点的数量以及每个拉格朗日点周围的欧拉点数量。在当前算法中，L在2D分析中为32，在3D分析中为33，但这取决于狄拉克δ函数(方程(4))中定义的外力的扩散率。相比之下，原始的边界条件强制IBM每次迭代需要N2次FLOPs。由于大多数情况下N?L，基于共轭梯度技术的IBM本质上比原始的边界条件强制IBM快得多。
基于共轭梯度技术的IBM的计算过程如下：该方法也尝试求解方程(16)中的Y，但方向是共轭的：
(21) Yk+1 = Yk+αkpk
(22) αk = (rk)TrkCk
(23) pk = rk?∑iO(10?)）时，每次迭代减少的FLOPs数量LN仍然会导致较大的计算成本。为了进一步应用于具有更高雷诺数或复杂几何形状的流动，需要在壁面周围获得更高的分辨率，因此需要更多的拉格朗日点。本文提出了一种新的策略，在壁面具有挤压几何形状时，避免显式评估A∈RN×N，从而进一步降低FLOPs。
在基于力校正的边界条件强制IBM的原始算法中，为了求解方程(7)、(8)中的校正算子，随着表示壁面的拉格朗日点数量N的增加，计算成本会按O(N2)的顺序增加。这是因为，当N增加时，矩阵A∈RN×N的大小也会增加，而迭代求解?的计算成本直接受到矩阵A大小的影响。如果能够减少矩阵A的大小，那么求解修正算子的计算成本也可以降低。因此，提出了一种新的策略来评估方程(8)，以减少计算成本。该算法的核心思想是利用这样一个事实：当几何形状在展向方向上是对称的时，方程(8)在对称位置会产生相同的解。这一概念在图2.1中进行了说明。由于对称位置的?解是相同的，因此可以避免重复计算。例如，?1和?5是对称的；因此，不必分别求解?1和?5，而是可以直接将?1的解赋值给?5。所提出算法的简化过程总结在图2.2中。当壁面包含挤压几何形状（如圆柱或等弦翼）时，可以通过简化矩阵A来降低计算成本。作为一个示例，考虑了一个在展向方向上周期性的等弦NACA翼型。图2.3显示了沿展向方向的壁面上拉格朗日点的分布，其中线条的交点对应于拉格朗日点的位置。所提出的算法适用于在展向方向上均匀分布或均匀拉伸的网格点，后一种情况的适用性将在后面讨论。在这项研究中，通过在每个垂直于展向方向的平面上重复相同的序列，从展向的一端到另一端来生成和排列拉格朗日点。

下载：下载高分辨率图片（376KB）
下载：下载全尺寸图片
图2.1. 所提算法的核心思想（q和r分别表示单个x-y平面上的拉格朗日点数量以及沿展向方向的对称平面数量）。对称位置的?解是相同的，因此可以在一个位置获得的解直接赋值给另一个位置。
下载：下载高分辨率图片（271KB）
下载：下载全尺寸图片
图2.2. 总体简化过程（N和q分别表示拉格朗日点的数量以及单个x-y平面上的拉格朗日点数量）。通过三步简化线性方程，可以减小矩阵大小，以提取必要的计算内容。
下载：下载高分辨率图片（412KB）
下载：下载全尺寸图片
图2.3. 所提算法适用的拉格朗日点配置示例（线条的交点表示拉格朗日点的位置）。理想情况下，几何形状在展向方向上应是对称的，并且沿展向方向的间隔应为均匀或均匀拉伸的，以获得最小的误差。
下载：下载高分辨率图片（639KB）
下载：下载全尺寸图片
图2.4. 考虑到?解的对称性时拉格朗日点分割的示意图。每种颜色代表同一对称平面上的拉格朗日点（q表示同一对称平面上的拉格朗日点数量）。每个平面上都有相同数量的拉格朗日点，q。

2.4.1. 简化方法1：利用?解的对称性
图2.4展示了表示翼型壁面的拉格朗日点的分布。这里，q表示单个x-y平面上的拉格朗日点数量，颜色（红色、蓝色和绿色）表示拉格朗日点位于同一平面上。这些点沿z方向均匀分布。由于拉格朗日点与欧拉网格之间的模式化对应关系，当第k个和第l个拉格朗日点具有相同的x和y坐标时，方程(7)中的每个Akl取相同的值。因此，这些拉格朗日点处的?值也必须相同。这种关系可以如下分析表达：(28)
?1,?2,?,?N=(?1,?,?q)r
这里，r是z方向上的欧拉或拉格朗日点数量，序列?N是通过将(r个?1,…,?q)块串联起来直到长度为N得到的。因此，虽然方程(7)、(8)中的?解原本有N个自由度，但通过考虑坐标中的模式化关系，自由度减少到了N/r (=q)。
(29)A11A12…A1NA21?????Aq1Aq2…AqN?1?2??N=11?1
这里，为了明确矩阵?中的元素数量为N，即使在[?1,…,?N]=[(?1,…,?q)]r的情况下，方程(29)中也使用了?N。

2.4.2. 简化方法2：利用矩阵A的稀疏性
图2.5显示了与求解任意x-y平面上?值相关的拉格朗日点。蓝色的拉格朗日点对于求解?的值是相关的。根据方程(4)，z坐标超出1.5Δz范围的蓝色拉格朗日点对应的Akl值为零。从这个角度来看，在图2.5中，只有与z方向上蓝色拉格朗日点相邻的红色拉格朗日点需要用来求解?的值。在分析层面，这可以通过以下过程表达。专注于垂直于展向方向的其中一个平面（例如pth平面）上的解，方程(29)现在等同于求解以下方程。(30)
0?0A1,(p?1)q+1?A1,(p+2)q
0…00?0A2,(p?1)q+1?A2,(p+2)q
0…0??????0?0Aq,(p?1)q+1?Aq,(p+2)q
0…0?1?2??N=11?1
这里，p∈{1,2,…,r?1}。这是因为在图2.5中灰色显示的(p?1)th之前和之后的平面（这些平面中的z坐标超过1.5Δz）与求解?无关，因此Akl可以被替换为0。由于简化后的矩阵在等式(30)左右两侧的对齐零值，它可以进一步转换成以下形式。(31)
0?0A1,(p?1)q+1?A1,(p+2)q
0…00?0A2,(p?1)q+1?A2,(p+2)q
0…0??????0?0Aq,(p?1)q+1?Aq,(p+2)q
0…0?(p?1)q+1??(p+2)q=11?1
下载：下载高分辨率图片（663KB）
下载：下载全尺寸图片
图2.5. 考虑到蓝色拉格朗日点解的对称性时相关拉格朗日点的示意图（蓝色：对其解感兴趣；红色：获取蓝色点解所需考虑的；灰色：对获取蓝色点解无关）。
然后，方程(29)可以修改为以下q×3q矩阵方程，并且在任何p(∈1,?,r?1)的情况下方程(32)都能得到满足。(32)
A1,(p?1)q+1A1,(p?1)q+2…A1,(p+2)q
A2,(p?1)q+1?????Aq,(p?1)q+1Aq2…Aq,(p+2)q
?(p?1)q+1?(p?1)q+2??(p+2)q=11?1

2.4.3. 简化方法3：用于Krylov子空间迭代方法的矩阵变换
Krylov子空间迭代方法能够在没有预处理的情况下用几次迭代就求解原始方程(8) [8]；因此，在这种情况下使用相同的迭代方法。上述简化将矩阵A的大小从原来的N×N减少到q×3q，但是A的非方形矩阵特性阻碍了Krylov子空间迭代方法的使用。为了应用Krylov子空间迭代方法，最好将矩阵A修改为方形矩阵，并利用以下关系：(33)
?(p?1)q+1,?(p?1)q+2,?,?pq=?pq+1,?,?(p+1)q
(34)
?(p+1)q+1,?(p+1)q+2,?,?(p+2)q=?pq+1,?,?(p+1)q
因此，(35)
?(p?1)q+1,?(p?1)q+2,?,?(p+2)q=(?pq+1,?,?(p+1)q)
3。
这里的矩阵右侧是通过将(?pq+1,…,?(p+1)q)块串联起来达到3q长度得到的。因此，方程(32)可以重写为(36)
A1,(p?1)q+1…A1,pq???Aq,(p?1)q+1…Aq,pq
?pq+1??(p+1)q+A1,pq+1…A1,(p+1)q???Aq,pq+1…Aq,(p+1)q
?pq+1??(p+1)q+A1,(p+1)q+1…A1,(p+2)q
?pq+1??(p+1)q=1?1
并且可以使用修改后的矩阵A来简化方程(8)。(37)
Amod∈Rq×q:
(37)Amod?mod=1
其中(38)
Amod=A1,(p?1)q+1+A1,pq+1+A1,(p+1)q+1…A1,pq+A1,(p+1)q+A1,(p+2)q
???Aq,(p?1)q+1+Aq,pq+1+Aq,(p+1)q+1…Aq,pq+Aq,(p+1)q
(39)?mod=?pq+1??(p+1)q
这些矩阵简化使得每次迭代只需要q2 (=N2/r2)次FLOPs，而不是原始边界条件强制IBM中的N2次FLOPs或共轭梯度技术基IBM中的LN FLOPs。例如，当展向方向的网格点数量为100时，每次迭代只需要原始FLOPs的1/1002。这种显著减少的FLOPs使得边界条件强制IBM能够应用于具有更多拉格朗日点的模拟中，特别是当壁面具有挤压几何形状时。组装修改后的矩阵Amod的实现过程总结在算法1中。计算方程(32)中的矩阵A的第一步与原始算法相同，但循环的大小减少了。方程(37)中获取矩阵Amod的第二步是新添加的。这里，p是从{1,…,r?1}中任意选择的，之后将?的解转移到其他平面上。由于组装矩阵Amod引起的计算开销与方程(37)中修正算子?的迭代收敛过程相比可以忽略不计（<0.01%）。

拉格朗日点也可以沿展向方向均匀拉伸以应用这些简化，从而扩大了应用范围。这是因为方程(7)中的矩阵A仅依赖于欧拉网格与拉格朗日点之间的相对距离。从分析上讲，这种依赖性可以通过以下方程明确表示。可以安全地假设欧拉网格和拉格朗日点的z坐标是对齐的，因为这是理想情况下定位拉格朗日点的最佳和最简单的方法。Pinelli等人[8]最初讨论了拉格朗日点与欧拉网格大小相同的理想距离。他们指出，将Δs/Δx减小到1以下会增加误差。因此，|zj?Zs|的解可以明确表示为以下模式之一。(40)
|zj?Zs|=Δzj0Δzj+q
这里，Δzj由|zj?zj?q|定义。由于每个垂直于展向方向的平面上包含q个拉格朗日点，Xj和Xj?q在z方向上是相邻的，并且具有相同的x和y坐标。当网格点沿展向方向以α的拉伸率均匀拉伸时，|zj?Zs|/Δzj的解可以写为(41)
|zj?Zs|Δzj=10α
因为Δzj+q=αΔzj。这表明矩阵A与展向方向的局部网格大小无关，因此可以像均匀网格情况一样进行简化。因此，即使网格在展向方向上被拉伸，方程(8)的解也不会改变。需要注意的是，在当前算法中，z方向上相邻的第二个拉格朗日点之间的拉伸率需要保持适度，以满足α2+α>1。尽管这在某种程度上限制了拉伸率（在方程(4)中，δh(r)考虑了外围的1.5Δz拉格朗日点），但可以通过α2+α>1来描述这一限制。即使拉伸率限制在82%到121.5%之间，这仍然在实际使用范围内是可以接受的。对于更高的拉伸率，可以通过考虑展向方向上的额外拉格朗日点以类似的方式修改矩阵A。例如，考虑展向方向上网格缩小侧的额外拉格朗日点，可以将拉伸率的限制从69%增加到145%，这是基于α3+α2+α>1得出的。对于更一般的拉伸率，考虑展向方向上的额外拉格朗日点不会导致明显的计算开销，因为额外的计算与矩阵A的大小无关。

此外，所提出的算法也适用于部分满足图2.3所示条件的壁面。当对称性部分被破坏时，只需要为那些拉格朗日点求解正确的?。例如，当等弦翼具有翼尖时，可以在翼面的展向内部应用所提算法，但在翼尖处不适用。因此，对于翼面的其余部分，应分别对那些更接近翼尖的拉格朗日点进行修正算子计算以获得?。简化和非简化区域之间的接口可以如下处理：在模拟开始时，根据坐标将每个拉格朗日点分类为简化区域或非简化区域。在每个时间步长中，首先获取简化区域中的修正算子，然后获取非简化区域旁边的修正算子的解，并将其作为ε的部分常数值。然后，可以使用原始算法来解决非简化区域中剩余的修正算子。在这里，可以在非简化区域计算修正算子时省略简化区域内修正算子的迭代计算，因为简化区域内的修正算子已经得到。这使得非简化区域内修正算子的额外计算成本也较小。3. 结果与讨论3.1. 计算效率与准确性分析为了评估边界条件增强型IBM的计算效率改进情况，通过求解不可压缩Navier-Stokes方程进行了折翼机的大涡模拟（LES）。计算域如图3.1所示。该翼型采用NACA0010配置，基于翼型的弦长C，计算域在x、y、z方向上的尺寸分别为12C、12C和2C，网格点数分别为1140、550和67。z方向上的展向长度为2C，边界条件在x、y方向上为出流条件，在z方向上为周期性边界条件。最小网格尺寸为(Δx,Δy,Δz) = (0.006C, 0.006C, 0.03C)。入口流在x方向的负端施加。雷诺数基于环境流速U和翼型的弦长，为3000。折翼机在y方向上具有正弦波振荡，振荡幅度为0.5C，频率为0.5U/C Hz，相对于环境流的攻角为零度。每个x-y平面上有332个（=q）拉格朗日点，整个域共有22,244个（=N）点。未解决的子网格尺度建模使用了动态Smagorinsky模型[18][19]。LES在Intel Xeon CPU Max 9480 @ 1.9 GHz上进行，并使用消息传递接口（MPI）进行并行化。图3.2显示了z方向上旋度的分布，该分布按翼型的振荡频率进行了归一化。下载：下载高分辨率图像（241KB）下载：下载全尺寸图像图3.1. 计算域（虚线箭头表示翼型的振荡运动）。图3.3展示了原始算法与改进算法之间的准确性和计算效率比较。图3.3(a)展示了修正算子?的比较，?按原始IBM算法的平均值进行了归一化。从翼型的不同位置采样了大约10,000个点，并绘制了时间步长。颜色刻度表示每个拉格朗日点获得的误差值（?simplified??original）/?originalˉ×100，其中?originalˉ是?original的平均值，最大差异幅度为0.058%，可以忽略不计。这里使用了稳定双共轭梯度（BiCGStab）线性求解器的迭代收敛策略[20]。值得注意的是，对于具有较强三维效应或不稳定流动或更高雷诺数的情况，修正算子的误差原则上也会保持类似的小幅度，因为修正算子仅取决于欧拉网格与拉格朗日点之间的相对距离，这与流动条件无关。下载：下载高分辨率图像（290KB）下载：下载全尺寸图像图3.2. 无量纲化旋度等值线。下载：下载高分辨率图像（280KB）下载：下载全尺寸图像(a)。准确性。下载：下载高分辨率图像（277KB）下载：下载全尺寸图像(b)。计算成本（对数刻度）图3.3. 比较(a)归一化的修正算子?（颜色刻度表示(?simplified??original)/?originalˉ×100），以及(b)原始算法与改进算法之间的计算时间（百分比表示IBM计算成本的占比）。“其他”表示排除IBM计算之外的LES计算成本。随后，比较了改进算法与原始算法的计算成本。比较时使用了880个和1760个CPU核心进行并行化。图3.3(b)以对数刻度显示了每个步骤的计算成本及其分解。这里，“其他”表示排除IBM计算之外的LES计算成本。使用原始算法的IBM的计算成本占总成本的94%至97%，而使用改进算法的IBM的计算成本仅占总成本的17%至24%，使得IBM计算的速度提高了52至160倍。实际上，计算速度的提高主要归因于上述修正算子计算的简化，这本身就使得速度提高了800至2800倍。因此，对修正算子?的评估和计算速度的提高表明，改进算法能够以显著高效和准确的方式再现?的解。值得注意的是，原始算法与简化算法之间在线性系统的条件处理上没有观察到差异。求解器容忍度保持不变，没有额外的预处理，迭代次数也没有明显差异。对于其他Krylov迭代方法，预计也会出现类似的计算效率和准确性效果，因为Krylov方法的计算成本主要来自涉及矩阵A的矩阵-向量乘法。预计减小矩阵A的尺寸也将降低任何Krylov迭代方法的计算成本。此外，图3.4展示了原始算法与改进算法之间的时间和时间平均升力及阻力系数的比较。时间和时间平均升力及阻力系数是通过积分翼型表面的压力和剪切应力得到的，这些数据是从欧拉网格插值到拉格朗日点获得的。原始算法与改进算法得到的时间和时间平均升力及阻力系数结果吻合良好，表明改进算法能够以显著高效的方式再现感兴趣的物理量，而不会损失准确性。下载：下载高分辨率图像（263KB）下载：下载全尺寸图像(a)。时间。下载：下载高分辨率图像（186KB）下载：下载全尺寸图像(b)。相位平均。图3.4. 原始算法与改进算法之间的(a)时间和(b)相位平均升力及阻力系数的比较。升力和阻力按动态压力进行了归一化，时间按翼型的振荡频率进行了归一化。3.2. 偏差敏感性分析改进算法基于拉格朗日点在展向方向上对称的假设，并在这种假设下求解修正算子。然而，在实际应用中，对称性假设可能在一定程度上被打破，需要研究偏离对称性假设对修正算子的敏感性。为了评估这种敏感性，使用原始算法和改进算法模拟了弦长变化的翼型周围的流动，并比较了结果。图3.5展示了用于敏感性分析的弦长变化翼型几何形状的示意图。该翼型采用NACA0010配置，展向跨度为60毫米，弦长从左端的基弦C线性变化到(a) 0.99C和(b) 0.90C，分别对应于图3.5中对称性假设偏离a=0.2度和a=1.2度。图3.6展示了原始算法与改进算法之间修正算子?的比较结果。颜色刻度表示每个拉格朗日点获得的误差值（?simplified??original）/?originalˉ×100，其中?originalˉ是?original的平均值，最大差异幅度分别为18.7%和81.0%，对应于对称性假设偏离0.2度和1.2度的情况。这表明改进算法对对称性假设的小偏差非常敏感。可能需要对修正算子的计算进行一些参数调整，以提高算法对抗对称性假设偏离的鲁棒性。下载：下载高分辨率图像（160KB）下载：下载全尺寸图像图3.5. 用于敏感性分析的弦长变化翼型几何形状的示意图。a表示基于展向方向的不对称度。下载：下载高分辨率图像（189KB）下载：下载全尺寸图像(a)。a=0.2度。下载：下载高分辨率图像（198KB）下载：下载全尺寸图像(b)。a=1.2度。图3.6. 原始算法与改进算法在不同对称性假设偏差下的修正算子误差。a表示图3.5中所示的基于展向方向的不对称度。下载：下载高分辨率图像（241KB）下载：下载全尺寸图像(a)。均匀缩放。下载：下载高分辨率图像（235KB）下载：下载全尺寸图像(b)。展向缩放。图3.7. 比较了(a)均匀增加分辨率或(b)增加展向长度l（粉红线和绿色虚线分别代表单独求解修正算子和整个IBM的计算成本，实心和空心符号分别代表原始算法和改进算法）。灰色条形表示排除IBM的LES计算成本。3.3. 可扩展性分析为了评估改进算法的可扩展性，将其计算复杂性与之前的算法进行了比较，如表1所示。这里，r和κ分别代表z方向上的欧拉点或拉格朗日点数量和条件数。边界条件增强型IBM对于挤压壁几何形状的复杂度为O(max(N,N2/r2))，这是因为原始修正算子的计算复杂度为O(N2)，由于矩阵A的大小从A∈RN×N减小到Amod∈Rq×q（q=N/r），因此简单减少了r×r。这意味着在改进算法中，如第2.4节所述的沿展向方向对称增加拉格朗日点数量，并不会增加修正算子的计算成本。然而，O(N)的基本计算复杂性无法避免，因为这种复杂性源于每个拉格朗日点的计算，例如从欧拉网格到拉格朗日点的插值或对外部假力的施加。通过使用单个处理器（Intel(R) Xeon(R) Silver 4210 CPU @ 2.20 GHz）在不进行并行化的情况下尝试模拟不同网格分辨率和展向长度下的圆柱体周围流动，研究了改进算法的定量可扩展性。在Pinelli等人[8]提出的原始算法中，欧拉点和拉格朗日点的分辨率理想情况下保持不变。因此，在此分析中，随着欧拉点数量的增加，拉格朗日点的数量也相应增加。图3.7展示了在不同(a)分辨率（均匀缩放）和(b)展向长度l（展向缩放）下，原始算法与改进算法之间的计算成本比较。均匀缩放以相同的速率增加了欧拉网格和拉格朗日点的分辨率。这里还明确显示了仅修正算子计算的计算成本。图3.7(a)中，当欧拉网格和拉格朗日点的分辨率增加时，改进算法中的修正算子计算成本始终比原始算法小103至104倍。此外，虽然在原始算法中修正算子计算占主导计算成本，但在改进算法中大部分计算成本来源于修正算子计算之外的计算，例如对每个拉格朗日点的速度场施加外部强迫项。这表明使用改进算法几乎消除了修正算子的计算成本，而这正是边界条件增强型IBM的缺点。请注意，在实际应用中，可以通过将欧拉网格和拉格朗日点分布在多个处理器上来进一步降低算法的计算成本，就像在原始算法中一样。表1. 边界条件增强型IBM算法之间的计算复杂性比较。IBM算法计算复杂性基于显式技术的IBM [15]O(N)多方向强迫IBM [6], [7]O(Nκ)基于共轭梯度技术的IBM [15]O(Nκ)边界条件增强型IBM [8]O(N2)用于挤压壁几何形状的边界条件增强型IBM（改进型）O(N),r≥N,ON2r2,r值得注意的是，该特性可以通过调整来提高横向方向的分辨率，同时保持恒定的横向长度。由于欧拉网格和拉格朗日点数量的增加，原始算法的计算成本会随着横向长度的增加而增加。另一方面，所提出算法的计算成本增加幅度没有原始算法那么大。实际上，修正算子计算的计算成本几乎保持不变。在(a)和(b)两种缩放分析中，简化后算法的整体计算成本与不包含IBM的LES的计算成本处于同一数量级。而原始算法的整体计算成本则比不包含IBM的LES高10到10000倍。这表明，IBM原本占据了整个模拟计算成本的主体，但简化后的算法显著降低了IBM的计算成本，使其与不包含IBM的LES处于同一数量级，这是边界条件强化IBM计算效率上的一个重大改进。不同IBM算法的可扩展性定量比较见表2所示。这些数值与Zhao等人[15]进行的圆柱体周围流动数值模拟结果进行了对比。所提出算法的计算成本在(a)均匀缩放和(b)横向缩放两种情况下进行了评估。其中，(b)横向缩放仅增加了横向方向的欧拉网格和拉格朗日点数量，从而提高了分辨率或域大小。所提出算法的计算成本在(a)均匀缩放情况下为O(N1.352)，在(b)横向缩放情况下为O(N1.008)。这两个数值都低于其他算法的计算成本，这些算法并未牺牲计算精度，例如多方向强迫IBM（O(N1.683)）、基于共轭梯度技术的IBM（O(N1.428）和原始边界条件强化IBM（O(N2.021)）。特别是(a)均匀缩放的结果表明，当横向方向的拉格朗日点数量增加时，计算成本没有呈指数级增长，与其他高精度算法不同。此外，所提出算法的(b)横向缩放结果与基于显式技术的算法（O(N1.014)）几乎相同，后者通过牺牲计算精度来提高计算效率。

**表2. IBM算法之间的计算成本缩放比较**

| 巧计算法 | 计算成本缩放 |
|---------------|-------------|
| 显式技术基础IBM | O(N1.014) [15] |
| 多方向强迫IBM | O(N1.683) [15] |
| 基于共轭梯度技术的IBM | O(N1.428) [15] |
| 边界条件强化IBM | O(N2.021) [15] |
| 边界条件强化IBM（针对挤压壁几何） | (a)均匀缩放 | O(N1.352) |
| | (b)横向缩放 | O(N1.008) |

**4. 结论**
提出了一种针对挤压壁几何的边界条件强化扩散界面IBM中修正算子计算的新简化公式。该公式采用了简化算法来求解修正算子，从而在保持计算精度的同时显著降低了计算成本。该算法利用墙体几何的对称性来简化耗时的矩阵计算，并对其计算效率和精度进行了评估。结果显示，整体IBM计算的效率提高了多达160倍，且解的误差不超过0.058%。此外，计算成本的缩放幅度始终低于之前提出的那些不牺牲计算精度的算法。

**cRediT作者贡献声明**
Manabu Saito：撰写原始草案、可视化制作、验证、方法论设计、数据分析与整理；
Ryoichi Kurose：撰写审查与编辑内容、提供监督、资源协调及概念构思。