自主开发CFD专用软件，流体力学方程组的建立，边界条件和初始条件的选择，以及收敛条件的分析都比较关键。基于对不可压缩流体流过孔板基本代数方程组的分析，本文自主开发了孔板流量计流场的数值模拟软件，分析了软件收敛的条件，并给出了层流和湍流条件下流出系数的计算结果，期望能为未来自主开发更大型的、针对特定产品的CFD专用软件起到促进作用。

1 数学模型

1.1 不可压缩流体轴对称流动基本方程组的无量纲化

以孔板流量计的管径d，平均流速um，d/um和ρu2m分别作为坐标、流速、时间和压力的特征物理量；x，r，u，υ，p，t的相对物理量分别为：xr=x/d，rr=r/d，ur=u/um，υr=υ/um，pr=p/ρu2m，tr=tum/d，可得到不可压缩流体轴对称流动基本方程组的无量纲形式

（1）

（2）

（3）

式中：Re=umd/ν；ν为运动黏性系数，ν=μ/ρ。

1.2 不可压缩流体轴对称流动基本方程组的收敛问题

国内许多CFD的研究人员之所以更乐意选用FLUENT等商业软件作为研究工具，除商业软件具有界面友好的优点之外，最重要的原因之一就是CFD基本方程组的收敛问题不好解决。在现有的CFD相关文献中，鲜见详细论述解决CFD基本方程组收敛问题的方法。基于此，本文给出了不可压缩流体轴对称流动基本方程组的收敛条件。

1.2.1 轴对称流动基本方程组的耦合问题

连续性方程（1）在形式上相当简单，但它与不可压缩流体的动量方程（2），（3）不耦合。为此，采取虚拟可压缩性法，即假定流体为可压缩的，但用不可压缩流的边界条件去求出收敛的解。收敛时，应有：?V→0。此方法的核心思想是将可压缩流体的连续方程化为

（4）

设中ρ为变量，而把式（4）中其他的ρ取为定值，即有

（5）

设虚拟可压缩流体p与ρ为最简单的线性关系，即

（6）

将式（6）代入式（5），可得

（7）

式（6）和式（7）中，为虚拟可压缩流体的音速。用时间相关法和动量方程（2）和方程（3）一起求解式（7），→0，就有??V→0，可获得不可压缩流体收敛的解。

1.2.2 Δt的选择

Δt的选择对于离散之后的方程组的收敛至关重要。选择既考虑对流（流动），又考虑扩散（黏性）的简单的Burgers方程作为分析的出发点，可得

（8）

式（8）中：ucr为可调的参数（ucr通常大于1），CRe=10n-2（Re=10n）为另一可调参数（假定Re在大范围变化时，Δtr基本上不变）。适当地增加ucr即可减少Δtr，以达到数值稳定性的要求。

1.2.3 虚拟可压缩性的选择

在式（7）中，虚拟音速a起着十分重要的作用，其取值适当与否，往往成为采用虚拟可压缩性法数值求解成败的重要因素，对此作了深入的研究。本文用以下方法选择

式中：β=d/D；d为孔板流量计的圆孔直径；D为孔板流量计的管径.实际的数值模拟说明，按式（9）取值时，可以得到收敛的结果。通常虚拟音速a远小于实际流体音速（如空气的实际音速在常温常压时约为330m/s）。

1.3 边界条件和初始条件的选择

基本方程组体现了一大类研究对象的共性，而边界条件和初始条件，则体现该研究对象的特性。共性与特性相结合，有可能给出该研究对象的解。因此，正确地选定边界条件和初始条件是数值模拟中十分重要的内容。

1.3.1 边界条件

边界条件为在求解域边界处（入口和出口截面、边界的固壁面处和对称轴线上）给定的p，u和υ之值。对于定常流动，这些边界条件是不变的。

1）入口和出口截面处u和υ的分布

对于不可压缩流体的层流流动，可设在入口和出口处，轴向流速u的分布为在圆管中充分发展层流的速度分布，即

（10）

出口处的相对轴向速度，也可按处理。对于轴对称流动，可设入口和出口处的υr=0。

2）在固壁面上u和υ的速度分布

对于黏性流动，显然有：

ur（管道固壁面）=0

υr（孔板两侧固壁面）=0

由于在本文中选择异点离散法，u和υ的求解点（离散点）均不在对称轴上，因此无需给出它们在轴线上的值。

3）边界上p值的给定

为体现定常流时p的特点，可给定入口处的p值，即p0=const。式中，p0为给定值，可设其为零。

1.3.2 初始条件

对于定常流，一般无需给定初始条件。但在本文所选择的时间相关法或时间推进法中，将定常流视为非定常流，并在与定常流相应的一些边界条件下求其收敛解（此解即为定常流的解）时，合适的初始条件则是不可少的。

1）p的初始条件

可设p0=0。

式中，上标“0”为t=0时刻。

2）υ的初始条件

可设υ0=0。

3）u的初始条件

按照进口和孔板孔径（大小为βD）处流量相等的条件，有

设孔板孔径处u为充分发展层流的速度分布，则有

设在孔板孔径中，u0的分布规律不变。对于u在流场中的分布，可分区按线性变化计算。

1.4 计算网格的划分策略

由于流经孔板的流动是轴对称的，只需在图1中的x-r平面上进行流场的数值模拟。由于流场固壁的边界形状相当简单，可采用沿x和r方向均布的网格分布，沿x方向，在孔板前有l1个网格，孔板范围内有i0个网格，孔板之后有l2个网格。改变l1值，可判断入口管道的合适长度；改变l2值，可判断出口管道的合适长度；改变i0值，可判断孔板厚度的影响。

图1 孔板流场网格的划分

在r方向，取m个网格。本文采用异点离散法，p，u和υ方程的离散点是不相同的。p，u和υ的（i，j）点位置如图2所示。

图2 p，u和v方程的离散点位置

2 计算结果

表1～表3为本文计算出的在层流条件下的流出系数（孔板厚度均选定为Δx）与试验值Ce的对比情况。

表1 β=0.4，m=20时试验值Ce与计算值Cc对比

表2 β=0.5，m=20时试验值Ce与计算值Cc对比

表3 β=0.6，m=20时试验值Ce与计算值Cc对比

孔板湍流流场的数值模拟是在其层流流场的数值模拟基础上，在动量方程中用（μ+μt）代替μ，同时在k－ε模式中，增加k和ε两个偏微分方程进行的。此时，连续方程不变，仅用时均值代替层流的变量值。表4、表5为本文计算出的湍流条件下的流出系数Cd与国际标准化组织（ISO）数据的对比情况。

表4 Cd计算结果与ISO数据的对比（0.2＜β≤0.5，Cd为角接取压时之值）

表5 Cd计算结果与ISO数据的对比（0.5＜β≤0.75，Cd为角接取压时之值）

3 结语

1）与采用FLUENT商业软件相比，自主开发特定产品的CFD计算软件具有设置灵活、针对性强的特点，而且可以更好地把积累的经验数据融合到计算软件中去，提高计算的准确性。为用商业软件FLUENT进行孔板流量计流场的数值模拟。显然，本文在计算精度上具有优势。

2）进行了层流条件下的孔板流量计流场的数值模拟，但是，只是对Re≤25的计算结果与试验数据的比较。本文则将层流条件下计算出的流出系数与给出的试验数据进行了比较系统的对比。不仅如此，与相比较，在同样自主开发软件的前提下，本文同时进行了层流条件下和湍流条件下的数值模拟。

3）国内许多CFD的研究人员之所以更乐意选用FLUENT商业软件作为研究工具，最重要的原因之一，就是CFD基本方程组的收敛问题不好解决。笔者期望本文总结出的这些经验能为未来自主开发更大型的、针对特定产品的CFD专用软件起些促进作用。同时，本文所述CFD软件采用FORTRON语言编写，功能比较单一，是针对孔板流量计流场数值模拟开发的，今后可以在提高软件界面的友好性以及功能性等方面作进一步改进。