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

Brianne Christopher 2018年 8月 20日

在建立流体流动仿真时,我们通常关注较大型系统中的单个(也可能是多个)组件,比如水处理厂中的泵或沉淀池。这自然而然地带来了一个问题:在不干扰过程的情况下,我们可以在多远的距离上应用边界条件?在本篇博客文章中,我们分析压缩率可忽略不计的均匀流体的内部和外部流动的入口和出口边界距离的效果。

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

CFD 仿真通常需要大量计算,我们自然会尽量减少仿真中的自由度。如果将自由度减少到极致,我们可能会得到一个入口边界和出口边界相交的几何结构。让我们设想一个横截面为半圆形的管道中的 90°管道弯头。

Geometry of a 90-degree bend in a pipe. 管道中 90 度弯头的几何结构。

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

如果使用上图显示的几何结构建立仿真,则入口和出口边界共享一条边。在许多情况下,单单这个现象就可能导致严重的收敛问题。不过,在这个特定示例中,经过几次迭代后,解会收敛。我们还考虑建立一个恰当的仿真,其中入口和出口管道延伸到十倍半径的长度(如下图所示)。

A 90-degree pipe bend with an extended inlet and outlet. 带延伸的入口和出口的 90 度管道弯头。

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

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

Side-by-side plots of the pressure variation across the pipe bend in COMSOL Multiphysics®. COMSOL Multiphysics® 中管道弯头上压力变化的并排绘图。
横截面为半圆形的管道中 90°弯头上的压力变化,其中绘制了上游边界上的压力表面图和管道壁上的等压线。左侧的绘图显示带延伸的入口和出口管道情况下的结果,右侧的绘图显示不带延伸管道情况下的结果。

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

(1)

c_{p}=\frac{\Delta p}{\frac{1}{2}\rho U^{2}},

不带延伸的入口和出口管道的情况下,值为 2.3,带延伸的入口和出口管道的情况下,值为 0.60。通过观察速度场可以了解更多信息。
An image of the velocity profiles and streamlines in the pipe bend. 管道弯头中的速度分布和流线图。

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

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

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

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

(2)

L_{E}=0.05D_{h}

\left(\frac{UD_{h}}

{\nu}\right)

适用于层流

(3)

L_{E}=4.4D_{h}\left(\frac{UD_{h}}{\nu}

\right)^

{1/6}

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

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

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

Schematic of the upper half of a plane. 平面上半部分示意图。
Schematic of a sharp 90-degree bend in a plane. 平面内 90 度急转弯头示意图。

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

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

(4)

u-iv=U/t,\hspace

{5mm}t=\sqrt{\frac{z-1}{z+1}},\hspace{5mm}

\zeta=i\log\left({\frac

{1+it}

{1-it}}\right)+\log\left({\frac

{1+t}

{1-t}}\right)

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

A plot of the pressure coefficient along the inner wall of the bend. 沿弯头内壁的压力系数图。

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

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

重力分析

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

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

An image showing the different phases in a sedimentation tank model. 显示沉淀池模型中不同阶段的图。

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

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

An image of a turbulent thermal plume model with inlet and outlet boundary conditions. 具有入口和出口边界条件的湍流热羽流模型图。

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

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

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

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

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

A plot of the turbulent flow around an airfoil. 机翼周围的湍流绘图。

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

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

A schematic showing the potential flow from an obstacle and its wake. 显示障碍物及其尾流的势流的示意图。

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

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

(5)

\begin{align}(u,\, v)&=(U+\frac{Qx}{2\pi R^{2}},\, \frac{Qy}{2\pi R^{2}}),\hspace{12mm}\text{in 2D} \end{align}

\begin
{align}(u,\, v,\, w)&=(U+\frac{Qx}{4\pi r^{3}}, \,\frac{Qy}{4\pi r^{3}}, \,\frac{Qz}{4\pi r^{3}}), \hspace{5mm}\text{in 3D} \end{align}

其中 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)

p+\frac{1}{2}

\rho\left|\bf

{u}\right|^2=p_\infty+\frac{1}{2}\rho U^2\Longrightarrow c_p=\frac{p-p_\infty}{\frac{1}{2}\rho U^2}=1-\left(\frac{\left|\bf{u}

\right|}

{U}

\right)^2

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

(7)

\begin{align} c_p&=-\frac{2}{\pi}\left(\frac{x}{R}\right)\frac{d}{R}+O\left(\left(\frac{d}{R}\right)^2\right) \hspace{5mm}\text{in 2D} \end{align}
\begin{align}
c_p&=-\frac{1}{4}\left(\frac{x}{r}\right)\frac{d^2}{r^2}+O\left(\left(\frac{d}{r}\right)^4\right) \hspace{6.2mm}\text{in 3D} \end{align}

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

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

Schematic of the potential flow around an obstacle with circulation in 2D. 二维模式下有环流的障碍物周围的势流示意图。
3D view of the potential flow around an obstacle with circulation. 有环流的障碍物周围的势流的三维视图。

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

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

(8)

\begin{align} (u,\,v)& = (U+\frac{\Gamma y}{2\pi R^2},\,-\frac{\Gamma x}{2\pi R^2}),\hspace{134.4mm}\text{in 2D} \end{align}
\begin{align} (u,\,v,\,w)& = (U+\frac{\Gamma y}{4\pi R^2}\left(\frac{z+s/2}{\sqrt{R^2+(z+s/2)^2}}-\frac{z-s/2}{\sqrt{R^2+(z-s/2)^2}}\right),\,-\frac{\Gamma x}{4\pi R^2}\left(\frac{z+s/2}{\sqrt{R^2+(z+s/2)^2}}-\frac{z-s/2}{\sqrt{R^2+(z-s/2)^2}}\right) \end{align}
\begin{align}
& -\frac{\Gamma (z+s/2)}{4\pi (y^2+(z+s/2)^2)}\left(1+\frac{x}{\sqrt{R^2+(z+s/2)^{2}}}\right)+\frac{\Gamma (z-s/2)}{4\pi (y^2+(z-s/2)^2)}\left(1+\frac{x}{\sqrt{R^2+(z-s/2)^{2}}}\right), \end{align}
\begin{align}
& \frac{\Gamma y}{4\pi (y^2+(z+s/2)^2)}\left(1+\frac{x}{\sqrt{R^2+(z+s/2)^{2}}}\right)-\frac{\Gamma y}{4\pi (y^2+(z-s/2)^2)}\left(1+\frac{x}{\sqrt{R^2+(z-s/2)^{2}}}\right)),\hspace{7mm}\text{in 3D}
\end{align}

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

(9)

\Gamma= \pi Uc(\alpha+\beta)

其中,\alpha是迎角,-\beta是“零升力”角(都以弧度为单位)。

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

(10)

\begin{align} c_p&=-\left(\frac{y}{R}\right)\frac{c}{R}(\alpha+\beta)+O\left(\left(\frac{c}{R}(\alpha\beta)\right)^2\right) \hspace{13.5mm}\text{in 2D} \end{align}
\begin{align}
c_p&=-\frac{1}{2}\left(\frac{y}{R}\right)\frac{s}{R}\frac{c}{R}(\alpha+\beta)+O\left(\left(\frac{s}{R}\frac{c}{R}(\alpha+\beta)\right)^2\right) \hspace{5mm}\text{in 3D} \end{align}

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

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

An image of a NACA airfoil simulation in 2D. 二维模式下的 NACA 机翼仿真图。

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

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

An image of the computational domain for a simulation of a stalling airplane. 用于失速飞机仿真的计算域图。

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

关于入口和出口边界条件定位的结论

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

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

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

后续操作

了解 COMSOL Multiphysics 附加模块“CFD 模块”中流体流动建模问题的相关特征和功能:

参考文献

  1. R.V. Churchill & J.W. Brown, Complex Variables and Applications, 5th ed., McGraw-Hill, 1990.

博客标签

CFD 模块 技术资料
正在加载评论...

博客分类


博客标签