CFD 仿真中如何设置入口和出口边界条件

在进行流体流动仿真时，我们通常只关注较大系统中的单个组件（也可能是多个），例如水处理厂中的泵或沉淀池。自然地，这引出一个问题：在不干扰过程的情况下，我们可以在多远的距离应用边界条件？这篇博客，我们将分析入口和出口边界的距离对压缩率可忽略不计的均匀流体的内部和外部流动的影响。

设置内部流动的入口和出口边界

CFD 仿真通常需要大量计算，为了减少计算成本并确保模拟的精度，我们会尽量简化几何模型和减少自由度。但如果过度简化，可能会导致不符合实际的情况，如入口边界和出口边界相交的几何结构。以一个横截面为半圆形的 90° 尖锐弯头管道为例来说明。在建模时，应确保入口和出口边界的设置合理，不相互交叉，且能准确反映流体在弯头内的流动特性。

横截面为半圆形的 90° 弯头管道。

如果模拟是基于上图显示的几何结构设置的，那么入口和出口边界共享一条边。在大多数情况下，仅这个问题就可能导致严重的收敛问题。不过，在这个特殊示例中，解在几次迭代后就收敛了。让我们再来考虑一个合理的模拟设置，其中管道的入口和出口延伸至长度为 10 倍管道半径的位置（如下图所示）。

入口和出口延伸的 90° 弯头管道。

基于水力直径 D_{h}=4A/P ，在雷诺数为 120 的情况下执行仿真，其中 A 是横截面面积， P 是横截面周长。在入口处应用均匀的速度分布，在出口处施加 0 法向应力。下图显示两个仿真中弯头处的压力，左侧是带延伸的入口和出口管道的情况，右侧是不带延伸管道的情况。对于带延伸的入口和出口管道的情况，从压力中减去下游边界上的平均压力，使两个结果在此边界上的平均压力都为 0。

横截面为半圆形的90° 弯头管道上的压力变化，其中绘制了上游边界上的压力表面图和管道壁上的等压线。左侧的绘图显示带延伸的入口和出口管道情况下的结果，右侧的绘图显示不带延伸管道情况下的结果。

仿真结果显示，不带延伸的入口和出口管道的情况下，压力变化要大得多。由于所应用的均匀速度分布和壁面无滑移 边界条件之间的不相容性，靠近入口的壁上存在急剧的压力梯度。左侧的绘图显示弯头上游侧的压力更加均匀，这表明当流动到达弯头时已充分发展。不过，压力并不是完全均匀的，在靠近锐角的位置，压力看起来略低，这表明弯头的上游效应。我们还看到上游边界对面的管壁上有一个滞流点。弯头的损耗系数定义为

(1)

不带延伸的入口和出口管道的情况下，值为 2.3，带延伸的入口和出口管道的情况下，值为 0.60。通过观察速度场可以了解更多信息。

横截面为半圆形的管道中 90° 弯头的速度分布和流线。

上图显示了弯头上游四个位置和弯头下游四个位置的速度分布以及中心平面的流线。在上游区域，我们可以看到均匀的速度分布如何演变成充分发展的速度分布。在弯头处，我们看到管壁上的滞流点朝向入口管道和相关的再循环区域。在急转弯头的下游还有一个再循环区域，我们可以看到，充分发展的速度分布首先在出口管道末端获得。所有这些在简单的几何结构（只包含 90° 弯头）中都不存在，我们得到错误的压降也就不足为奇了。

入口 和出口 边界特征中提供的充分发展的流动 选项，可用于避免入口和出口管道过长。前面两个图中的结果充分表明，我们应该在距离弯头一定距离处应用这些条件以获得好的结果。那么我们需要在上游和下游多远的位置应用充分发展的流动 选项？如果我们将入口管道和出口管道分别从弯头处延伸 1 个半径的长度，弯头的损耗系数将变为 0.54，而如果每个方向延伸 2 个半径的长度，则损耗系数变为 0.58。从此时开始，收敛到 0.60 的速度会变慢。因此，在本例中，每个方向延伸 2 个半径的长度是不错的权衡方法。

随着雷诺数的增大，弯头下游的再循环区域长度将增大，最终变得不稳定。对于雷诺数为 1200 的情况，如果管道末端应用充分发展的流动 选项，当出口管道延伸超过 20 个半径长度时，损耗系数不会有明显变化。根据管道入口长度的相关性，

(2)

适用于层流

(3)

适用于湍流，我们可以估计在弯头下游多远处可获得充分发展的流动剖面。请注意，湍流入口长度通常比高雷诺数层流入口长度短。入口长度必须达到雷诺数 O(10^8)，才能达到水力直径 O(10^2)。

对于雷诺数为 120 和 1200 的两种层流情况，从 (2) 获得的入口长度分别约为半径的 7.5 倍和 75 倍。通过在出口处使用充分发展的流动 选项，我们获得了良好的结果，出口管道相当于这些长度的 1/3。

上游效应将随着雷诺数的增大而减小，这是因为，纳维-斯托克斯方程的椭圆特性会随着雷诺数的增大而减弱。我们可以通过观察类似几何结构的势流来估计上游效应区域。

使用 Schwarz-Christoffel 转换法将上半个平面映射到 90° 急转弯头。

使用 Schwarz-Christoffel 转换法（参考文献 1），复z平面的上半个平面可以被映射到复 \zeta 平面的 90° 急转弯头。入口位于 \zeta 平面的 -i\infty 处，对应于 z 平面原点的源，而两个平面的出口都位于 \infty 处。\zeta 平面中弯头的外角和内角分别对应于 z 平面中的点 -1 和 1。\zeta 平面内的速度场以隐式形式获得

(4)

下图显示沿内壁的压力系数随势流解的弯头上游无量纲距离变化的情况。

沿 90° 急转弯头内壁上游的压力系数。

图中，压力系数基于局部压力与远上游压力之差，h 是通道宽度。我们发现，当上游到弯头的距离为两个管道宽度时，压力系数为 O(10^{-3})。因此，如果我们在入口处使用充分发展的流动 选项，我们只需将入口管道（或通道）向上游延伸两三个水力直径的长度。

重力分析

在模型中包含重力的情况下，入口 和出口 边界特征中的充分发展的流动 选项带有附加的静水压力补偿 选项（不可压缩流动）或静水压力补偿近似 选项（弱可压缩或可压缩流动）。附加选项给出了不可压缩流动边界上的精确静水压力分布，以及弱可压缩流动和可压缩流动的良好近似。当入口或出口边界处的流动明显分层时（例如在多相流中），必须格外小心。在这些情况下，可能有必要增加一个腔室，使流动平行于重力矢量。

在对抗重力时也可能出现问题。下图显示了一个大型沉淀池，其停留时间较长，其中允许悬浮（重）相沉淀并通过底部出口流出，轻相通过靠近外缘的环形出口垂直排出。灰色流线对应于轻相的速度场，而黑色流线对应于重相的速度场。重相的一小部分通过轻相出口排出。这里，重相沿与重力相反的方向流动，当一些悬浮颗粒再次下落时，在外缘附近形成一个小涡流。这个小涡流会对时间步产生负面影响，导致总计算时间较长。一个可能的补救办法是添加溢流孔（溢流堰），使流动沿重力方向流出。

沉淀池中分散相的体积分数（彩色图）和流线，灰色表示轻相，黑色表示重相。

另一个流入和流出边界相交的示例发生在模拟热羽流时，如下图所示。此例中，在流入边界（图中的圆柱面）没有指定入口 边界条件，而是应用了开放边界 特征。在开放边界 特征和出口 特征（顶部边界）中，都使用了静水压力补偿近似 选项。这是必需的，原因是，模型中浮力引起的压力比静水压力小三个数量级。另一个同样重要的选项是出口 特征中的抑制回流 选项。

湍流热羽流，显示静水条件下压力偏差的速度大小（彩色图）和等值线。

顶部边界等值线的微小扰动是由于与恒压的不一致性造成的，这些扰动可以通过在非等温流动 多物理场耦合节点中使用布辛涅斯克近似 选项来消除。

放置外部流动的入口和出口边界

在外部流动应用中，例如车辆和建筑物周围的流动，远离障碍物的条件通常设置为入口边界上的恒定速度矢量和出口边界上的恒定压力。那么问题又来了，与应用这些条件的障碍物的距离会在多大程度上影响解？对于外部流动，事实证明该距离随模型的空间尺寸而变化。对于二维模型，所需距离比三维和二维轴对称模型大一个数量级。我们再次研究理想的势流解，尝试理解其中的原因。

障碍物周围的外部流动在固体表面的边界层中产生涡流。障碍物不同侧面上的边界层可能在后缘汇合，形成一个在下游用平流输送成尾流的薄涡旋面。如果任何一侧的边界层由于不稳定性或尖锐凸角的存而与障碍物分离，尾流将会更宽。在任何一种情况下，流向下游的涡流都被限制在尾流内，尾流外的流动近似于无旋流动。

NACA 机翼周围的湍流。剖面上侧的边界层在后缘之前分离。

障碍物及其尾流取代了自由流的流线，我们可以将离障碍物很远的流动看作是均匀流动与源的总和。

远离障碍物的势流及其尾流。

二维和三维模式下产生的速度场可以表示为

(5)

其中 R=\sqrt{x^{2}+y^{2}}, r=\sqrt{x^2+y^2+z^2}，源位于原点，自由流则流向正 x 方向。

在任何一种情况下，源强度都与障碍物的大小有关。在源的 x 位置，流线的位移在二维模式下是 y_{0}=Q/(4U)，在三维模式下是 r_{0}=\sqrt{Q/(2\pi U)}。

下游的极限值为别是 y_{\infty}=Q/(2U) 和 r_\infty=\sqrt{Q/(\pi U)}。出于当前的估计目的，我们可以使用二维模式下的 2y_{0} 和 2y_\infty 以及三维模式下的 2r_{0} 和 2r_\infty 作为障碍物大小的代表值。其中，我们在二维模式下使用 d=Q/(2U)，在三维模式下使用 d=\sqrt{2Q/(\pi U)}。根据伯努利方程，我们得到了远距离处压力系数的估算值

(6)

将势流速度场与源强度的估算值一起代入，我们可以得出

(7)

因此，压力系数在二维模式下减小为 d/R，在三维模式下减小为 (d/r)^{2}。为了减少外部边界条件（比如 O(10^{-2})）的影响，在二维模式下，我们必须将计算域的外部边界定位在 100 个障碍物大小的距离之外，在三维模式下，定位在 10 个障碍物大小的距离之外。

根据障碍物的形状和方向，涡旋脱落可能导致产生侧向力（升力）的环流。离障碍物很远的势流可以用均匀流和点涡（二维模式下）或均匀流和马蹄形线涡（三维模式下）来近似。

有环流（升力）的障碍物周围的势流。在二维模式下（左），势流包含 x 方向的均匀流和位于原点的点涡。在三维模式下（右），势流包含 x 方向的均匀流和马蹄形涡流，马蹄形涡流在 z 方向跨度为 s，在 x 方向延伸到无限远处。

在距离障碍物很远的位置，对应于上图的势流速度场由下式给出

(8)

请注意，通过将 z 设置为零并让 s\rightarrow\infty，我们可以根据三维解获得二维解。在大多数可实现的情况下，环流可以通过下式与障碍物的自由流速度和流向尺寸（弦）c 相关联

(9)

其中，\alpha 是迎角，-\beta 是“零升力”角（都以弧度为单位）。

后者源于障碍物的形状（曲率）；比如，机翼的弧度。将渐进势流解和 \Gamma 的表达式插入压力系数的定义中，得到

(10)

总偏转角 \alpha+\beta 必须至少比 1 个单位小一个数量级，这样 \Gamma 的估算值才成立。对于球体，维度 d、c 和 s 相等。因此，环流对外部边界附近的限制不如源的限制严格。

对于机翼来说，三个维度的数量级各不相同 s\sim 10c\sim 100d。在二维模式下，由于 d\sim c(\alpha+\beta)，点涡造成的限制与源产生的限制大小相等。如果要对三维机翼进行单独建模，那么线涡造成的限制将比源造成的限制严重 100 倍。通常，机翼通过 d\sim c 附在机身上，在这种情况下，三维模式下两种限制的数量级相同。下图显示 14° 迎角下 NACA 0012 机翼周围流动的二维仿真。为了将外部边界条件的影响最小化，我们在每个方向上将域延伸 100 个弦长。在该示例中，相关长度比例为 c\alpha，对于对称剖面，\beta=0。根据以上估算，压力系数为十分之几。

14° 迎角下 NACA 0012 机翼的二维仿真。

下图显示失速飞机在 20° 迎角下的三维仿真。计算域由半径为 15 米的半球和高度为 30 米的圆柱体界定。翼展约为 18 米，机身直径为 2.4 米，最大弦长乘以迎角约为 1.3 米。插入这些数值得出的压力系数为百分之几，该值偏高。因此，如果域延伸到离飞机更远的地方，仿真结果可能会有所改善。

根据速度大小着色的计算域，用于失速飞机的仿真。

关于设置入口和出口边界条件的总结

本文我们介绍了使用理想流动理论和经验关系式来确定流入和流出边界的合适位置。对于内部流动，我们使用层流和湍流的经验关系式来确定获得充分发展的流动所需的管道长度，相应地向上游和下游扩展域能够显著提升流动仿真的精度。不过，这个长度随着雷诺数增大而增加，尤其是对于高雷诺数层流而言，长度会变得过长。在入口 和出口 边界特征中使用充分发展的流动 选项可大大缩短这一长度。在下游，长度减小三倍似乎可以合理平衡精度和计算成本。使用势流理论表明，上游距离不需要超过多倍水力直径的长度。充分发展的流动 选项在重力处于活动状态时也会分析重力，但是当流动明显分层时（例如在多相流中），仍可能出现问题，在这种情况下，建议将出口转向重力方向。

对于外部流动，势流理论用于估计在多远的距离范围内可以忽略不计障碍物周围流动引起的压力变化。我们发现，在二维模式下，压力根据 \sim D/R 而变化，在三维和二维轴对称模式下，压力根据 \sim A/R^2 而变化，其中，R 是到障碍物中心的距离，D 和 A 分别是投射到与自由流正交的平面上的障碍物的长度和面积。

希望这些评估结果在您建立自己的仿真时有所帮助，不过还是要记住验证您的结果。当对内部流动使用充分发展的流动 选项时，最简单的方法是改变入口和出口通道或管道的长度，查看结果是否发生变化。对于外部流动，查看速度是否偏离压力边界上自由流速度场所容许的公差（尾流除外），对速度边界上的压力也执行类似操作。

编者注：这篇博客更新于 2023/4/25。用于流体仿真的 充分发展的流动边界条件已经是 COMSOL Multiphysics^® 核心模块的一个可选功能了。