第2章　放射性核素输运的数学模型

2.1　流动模型

2.1.1　自由面跟踪控制方程

为了模拟河道中的水流运动，首先需要考虑的问题就是如何跟踪与捕捉作为动边界的自由水面的位置。本书采用在海洋工程研究领域已大量使用的VOF方法（Hirt&Nichols，1981）来模拟河道中水面位置的变化。标准的VOF方法采用如式（2.1）的控制方程来分辨水面位置（Hemida，2008）。

式中：α为体积分数；U为流速。

α=1代表该网格单元为水体单元，α=0代表该网格单元为气体单元，0＜α＜1代表网格单元为水、气混合单元。在VOF方法中，单元密度ρ与动力黏度μ的表达见式（2.2）和式（2.3）。

式中：ρw为水体密度；ρa为气体密度；μw为水的动力黏度；μa为空气的动力黏度。

为了提高VOF方法对自由表面的捕捉精度，近来出现了一种改进的VOF方法。该方法在标准VOF方法基础上，引入压缩项来提高界面捕捉精度。改进后的VOF方法中，体积分数α满足式（2.4）的控制方程（Hemida，2008）：

式中：Ur为压缩速度。

数值计算中，压缩项，即式（2.4）左边第三项只需在界面区域进行处理。本书采用改进的VOF方法来跟踪河道水面位置的时空变化。

2.1.2　三维流动控制方程

对于弯曲较大、断面突变或者丁坝较多的天然河道，水流运动十分复杂，流线的不平行程度或者弯曲程度变化剧烈，垂直于流线的方向会产生较大的惯性力，静压假定不成立，准三维模型失效，需建立真三维模型来描述。

对于放射性核素在河道中的输运模拟而言，可考虑河道中的水流以及河道上方的空气为不可压缩流体。于是，控制水体和空气运动的连续性方程和动量方程见式（2.5）和式（2.6）（Versteeg&Malalasekera，2007）。

式中：p为压力；μ为动力黏性系数。

当α=1时，ρ=ρw，μ=μw，式（2.5）和式（2.6）就变为水流运动的控制方程。当α=0时，ρ=ρa，μ=μa，式（2.5）和式（2.6）就变为河道上方空气流动的控制方程。

真实河道中的流动在绝大多数情形下为湍流，包含了大量不同尺度的涡结构，且涡强度与方向都是随机的。湍流的高度非线性特性使得通过直接求解式（2.5）和式（2.6）来模拟湍流的全部细节变得极其困难。实际上，湍流的这些细节对于解决实际工程问题意义并不太大。对于河道中放射性核素输运的预报、诊断与辨识，更多关注的是湍动引起的流场的平均效果。本研究将采用雷诺平均法来研究河道中的湍流。

根据雷诺平均法的思想，任意时刻的瞬时速度U、压力p可进行的分解见式（2.7）和式（2.8）（Versteeg&Malalasekera，2007）。

于是式（2.5）和式（2.6）可用时均速度和时均压力表达，见式（2.9）和式（2.10）

为了表述方便，后续所有时均物理量的上画线均省略。于是可得

雷诺平均后出现了雷诺应力项，使得式（2.11）和式（2.12）并不封闭，需采用湍流模式进行模化。

2.1.3　湍流模型

目前，绝大多数放射性核素的输运模拟均采用涡黏模型。Margvelashvili（2000）采用零方程湍流模型模拟了切尔诺贝利核电站冷却池和第聂伯河口中的水流运动。该模型通过引入Richardson数来反映流动分层对涡黏系数的影响，并取临界Richardson数为10。零方程模型属局部平衡假设模型，无法考虑湍动的对流和扩散效应，对回流与分离流不适用。同时，复杂流情形下的混合长度难以确定，影响其工程应用的推广。

采用一方程模型模拟了与核素输运相关的湍动羽流。此外，在THREETOX核素输运模型中，垂向涡黏系数也通过一方程模型来获取（，2005a）。一方程模型的涡黏系数采用湍动能和混合长度来表示，克服了零方程模型局部假设的局限性。然而，一方程模型中的混合长度仍然需通过经验关系来确定，很难找到一个普遍有效的计算公式。

在两方程湍流模型中，涡黏系数由湍动能和湍动耗散率描述，解决了零方程模型与一方程模型中混合长度难以确定的问题，具有较强的适应性，在水动力学、气动力学研究领域已得到广泛应用。两方程模型中的标准k-ε模型，可对大量工程关注范围尺度上的涡流结构进行准确分辨，是目前航空、航天、海洋工程等诸多研究领域中应用最为广泛的两方程湍流模型。在现代湍流模拟中，虽然还存在大涡模拟、谱方法等其他湍流模拟方法，但其要求的计算量非常大，对实际河道中放射性核素输运模拟并不合适。

本书采用标准k-ε模型来模拟雷诺应力项。该模型基于Boussinesq假定，将雷诺应力与时均速度梯度相关联（Versteeg&Malalasekera，2007）：

式中：μt为湍动黏度；k为湍动能；I为单位二阶张量。

湍动黏度是湍动能和湍动耗散率ε的函数：

式中：Cμ为经验常数，取Cμ=0.09。

湍动能k与湍动耗散率ε由控制方程所描述（Versteeg&Malalasekera，2007），见式（2.15）和式（2.16）：

式中：νkeff为湍动能的有效扩散系数；νεeff为湍动耗散率的有效扩散系数；G为时均速度梯度引起的湍动能产生项，C1ε=1.44，C2ε=1.92。

νεeff和G可分别表达为

于是可得

式中：μeff为有效黏度，可由式（2.21）计算：

为了计算方便，引入具有压力量纲的物理量prgh，即

于是，式（2.20）可变为（Rusche，2002；Maki，2011）：

2.1.4　河道近壁区流动模型

大量试验表明，对于有壁面作用的充分发展湍流，近壁区流动可分为三个子层：即黏性底层、过渡层、对数率层（White，2006）。标准k-ε模型为高雷诺数模型，只适用于充分发展的湍流核心区，对于河床等近壁区域，并不适用，需进行特殊处理。对于近壁区流动，目前有两种处理方法：一是在近壁区使用低雷诺数k-ε模型；二是在近壁区使用壁面函数法。第一种方法需要将壁面附近网格进行局部加密，以分辨黏性底层与过渡层的流速分布，计算量大。第二种方法，即壁面函数法，只需保证近壁面网格的体心位于对数率层即可，计算量较第一种方法小。本书采用壁面函数法来处理河床近壁区流动。

真实河床底部并非光滑，底部粗糙会对流动产生影响。考虑粗糙度影响的河床底部湍流运动黏度采用式（2.24）计算（苏军伟，2012）：

式中：ν为运动黏性系数；κ为卡门常数，κ=0.41；E为经验常数，E=9.8；ΔB为反映床面粗糙对湍动黏性的影响。

在Open FOAM开源代码中，y*的分界值取10.97，而在Fluent软件中，y*分界值取11.225。y*定义见式（2.25）：

式中：kw为离壁面最近网格单元中心的湍动能。

ΔB与河床底部沙粒的形状和大小有关，可采用式（2.26）计算（Cebeci&Bradshaw，1977）：

式中：Cs为粗糙度常数，取值范围为[0，1]；为无因次壁面粗糙高度，定义为

式中：Ks为粗糙高度。当时，黏性底层较厚，粗糙高度完全淹没在黏性底层中，管壁粗糙对湍流结构无影响，即水力光滑区。当时，黏性底层较薄，粗糙高度部分深入到湍流核中，影响湍流结构，即湍流过渡区。当时，黏性底层已变得很薄，管壁粗糙高度完全进入湍流核中，即水力粗糙区。