多相流数值模拟常用数学模型

第6章　多相流数值模拟常用数学模型

解决多相流问题的第一步，就是挑选出最能符合所要解决的实际流动情况的多相流模型。确定各相的自身特性以及相与相之间（包括气泡、液滴和粒子）耦合的程度，针对不同的使用背景选择恰当的数学模型。本章将重点介绍目前多相流数值模拟所常用的数学模型及其使用的背景，而运用这些模型进行数值计算的案例将在第8章进行介绍。

6.1　无滑移单流体模型

6.1.1　模型的应用背景及特点

无滑移单流体模型的构建基于欧拉框架。在气固两相流动系统中，如果假定颗粒速度与当地气体速度相同，颗粒的湍流扩散系数与当地气体的湍流扩散系数也相同，且颗粒温度与流体温度相同或等于常量，则可将气固两相合并为一特殊的单相流体，其控制方程可以大大简化。

这种模型由于形式简单，计算方便，在早期两相流研究中使用较广。由于在大多数情况下悬浮体中的流体与颗粒之间总会存在速度差和温度差，特别是在起始段差别可以很大，因此应用这种方法算出的结果与实际差别较大。只有当颗粒尺寸很小时，才能得到满意的结果。以工厂实际问题利用单流体模型对粒度很小的钛白粉（dp＜0.5μm）的沉积计算为例，其结果与实际测定值比较接近。

6.1.2　主要数学模型

或

以及

无滑移单流体模型目前已较为成熟，当采用该模型进行数值计算时，可以自行编程或借助商业计算流体力学软件，例如FLUENT，CFX等。

6.2　双流体模型

6.2.1　模型的应用背景及特点

双流体模型是基于欧拉框架进行构建的，是一种相对较为复杂的多相流模型，它将各相看成相互渗透、耦合但又具有不同运动特征的连续介质，建立了一套包含有两组动量方程和连续性方程的方程组来求解每一相。压力项和各界面交换系数是耦合在一起的，耦合的方式则依赖于所含相的情况，颗粒流（流－固）的处理与非颗粒流（流－流）的处理是不同的。对于颗粒流，可应用分子运动理论来求得流动特性。不同相之间的动量交换也依赖于混合物的类别，该模型的应用包括气泡柱、上浮、颗粒悬浮以及流化床。比之单流体模型，双流体模型考虑了固相的湍流输运以及气固两相间相互滑移引起的阻力，其计算结果在大多数情况下更接近于实际。该模型可以对各相进行单独的计算，每一相都有单独的守恒方程，具有很大的适应性。但由于方程数的增加使计算趋于复杂，特别是当颗粒结构尺寸存在差别时，各相都要进行独自计算迭代，计算量是巨大的。此外颗粒相的扩散系数、导热系数和粘度的确定也仍然存在问题，并且双流体模型在计算中还可能产生伪扩散。

6.2.2　主要数学模型

6.2.2.1　守恒方程

双流体模型的守恒方程也是由质量守恒方程（连续性方程）、动量守恒方程、能量守恒方程组成，详细描述如下：

1）连续性方程

q相的连续性方程为

这里，α为q相的体积分数，ρ是q相的物理密度是q相的速度表示了从第p相到q qq相的质量传递。从质量守恒方程可得

和

2）动量守恒方程

其中，τq是第q相的压力应变张量（stress－strain　tensor）

方程（6－9）中，μq、λq是q相的剪切和体积粘度是外部体积力是升力是虚拟质量力，其具体数学表达式可参见本书第4章是相之间的相互作用力；p是所有相共享的压力。

是相间的速度，定义如下：如果（也就是，相p的质量传递到相q；如果m　·pq＜0（也就是，相q的质量传递到相p）

方程（6－9）必须有合适的表达使相间作用力R→pq封闭。这个力依赖于摩擦力、压力、内聚力和其他影响，并服从条件

可使用下面形式的相互作用力：

这里，Kpq（＝Kqp）是相间动量交换系数。液－液两相间动量交换系数可表示为Kll，对于颗粒流动，液－固交换系数为Ksl，固－固交换系数为Kss。下面将分别进行介绍。

6.2.2.2　相间交换系数

1）液－液交换系数

对于液－液流动，每个第二相被假定为液滴或气泡的形式。如何把流体中的一相指定为颗粒相有着重要的影响。例如，流动中有不同数量的两种流体，起支配作用的流体应作为主要流体，因为稀少的流体更可能形成液滴或气泡。对于这些气泡，液－液或气－液混合类型的交换系数可以写成以下通用形式：

这里，颗粒松弛时间τp定义为：

式中，dp是p相液滴或气泡的直径。

曳力函数f对不同的交换系数模型而言定义不同，但几乎都包含一个基于相对雷诺数（Re）的曳力系数（CD）。

（1）Schiller　and　Naumann模型

这里

Re是相对雷诺数。主相q和第二相p的相对雷诺数从下式获得

第二相p和r的相对雷诺数从下式获得

这里，μrp＝αpμp＋αrμr是相p和r的混合速度。

（2）Morsi and Alexander模型

这里

Re由方程（6－16）和（6－17）定义，a1，a2，a3定义如下：

Morsi　and　Alexander模型是最完善的，它频繁地在雷诺数的大范围内调整函数定义，但是这个模型使用起来比其他模型更不稳定。

（3）对称模型

这里

对于此模型，在流动中，区域内的某个地方的第二相（分散相）在另一个区域变成主相（连续相）。例如，如果空气注入充满一半水的容器的底部，在容器的底半部空气是分散相，在容器的顶半部空气是连续相。这个模型也用于两相之间的相互作用。

2）液－固交换系数

液－固交换系数Ksl以下面的通用形式写出：

颗粒的松弛时间τs定义为

这里，ds是（s）相颗粒的直径。

曳力函数f在不同的交换系数模型中其定义是不同的。

（1）Syamlal－O’Brien模型

这里，曳力函数采用由Dalla　Valle给出的形式

这个模型基于流化床或沉淀床颗粒的末端速度的测量，并使用了体积分数和相对雷诺数的函数关系式

这里，下标l是液体相，s是固体相，ds是固体相颗粒的直径。

此时，液－固交换系数有如下形式：

这里，vr，s是与固相相关的末端速度

当固相的剪切应力根据Syamlal等人提出的定义时［参见方程（6－46）］，这个模型是合适的。

（2）Wen－Yu模型

在该模型下，液－固交换系数有如下形式：

式中，CD由下式确定：

该模型适合于稀相系统。

（3）Gidaspow模型

该模型是Wen－Yu模型以及厄冈（Ergun）方程的联合。

当αl＞0.8时，液－固交换系数Ksl有如下形式：

式中CD为

当αl≤0.8时

对于密集的流化床，建议使用这个模型。

3）固－固交换系数

固－固交换系数Kss有如下形式：

式中，ess为恢复系数；Cfr，ss为第p和第q相之间的摩擦系数；dp为固体p颗粒的直径；g0，ss为径向分布函数。

6.2.2.3　固体压力

对于可压缩机制下的颗粒流动（也就是，固体的体积分数小于允许的最大值的地方），固体压力应独立计算并且用做颗粒相动量方程中的压力梯度相▽ps［参见式（6－9）］。因为麦克斯韦（Maxwellian）速度分布用于颗粒，将颗粒拟温度引入了模型并出现在固体压力和粘度的表达式中。由于颗粒的碰撞，固体压力由动能项和第二相组成

式中，Θs为颗粒拟温度，该值与颗粒运动的波动动能是成比例的，将在本部分的后面描述；函数g0，ss控制了从α＜αs，max（这里固体颗粒之间的距离可以继续减小）的可压缩条件到α＝αs，max（这里距离的进一步减小不会发生）的不可压缩条件。

径向分布函数g0是一个当固体颗粒相变密时用于修改颗粒之间碰撞概率的修正因子。这个函数也可解释为小球之间的无量纲距离

式中，s是颗粒之间的距离。从方程（6－42）可以观察出对于稀疏固体相，s→∞，所以g0→1。当固体相紧凑到一定限度内，s→0，g0→∞。径向分布函数与关于非均匀气体由Chapman和Cowling提出的理论的χ因子紧密联系。对于稀有气体，χ等于1，当分子靠得非常近以致运动不可能发生时，它会逐渐增加并趋向无穷大。

径向分布函数可由下式计算：

当固体相数大于1时，方程（6－43）扩展为

这里，αl，max是由你在问题的设置过程中指定的，并且

6.2.2.4　固体剪切应力

固体应力张量包含由于平移和碰撞而从颗粒的动量交换中产生的剪切和体积粘度。粘度的摩擦分量也可以包含在当固体颗粒相达到最大固体颗粒分数时出现的粘塑性变迁中。

固体剪切粘度包括碰撞部分、动力部分和摩擦部分

1）碰撞粘度

剪切粘度的碰撞部分为

2）动力粘度

Syamlal等人提出的表达式

Gidaspow等人提出的表达式

3）体积粘度

固体体积粘度解释为颗粒压缩和扩张的抵抗力，根据Lun等人的研究，有以下形式：

4）摩擦粘度

在低剪切密集流动中，固体的第二相体积分数接近于压缩极限，应力的产生主要是由于颗粒之间的摩擦。

如果计算中包含摩擦粘度，可使用Schaeffer的表达式

式中，ps是固体压力，φ是内部摩擦角，I2D是偏应力张量的第二不变式。它也可以被指定为常数或用户定义的摩擦粘度。

6.2.2.5　颗粒拟温度

第s固体相的颗粒拟温度是与颗粒的随机运动的动能成比例的。从动能理论得到的输运方程采用如下形式：

式中，为固体应力张量产生的能量；kΘs▽Θs为颗粒能量扩散通量（kΘs是扩散系数）；γΘs为能量的碰撞耗散；φls为第l液体相或固体相和第s固体相之间的能量交换。

能量的碰撞耗散γΘs代表了由于颗粒之间的碰撞在第s固体相内的能量耗散率。这项也可以由Lun等人提出的表达式描述

从第s固体相到第l液体相或固体相粒子速度的随机波动动能的传递由φls描述

6.2.2.6　湍流模型

为了描述单相中速度及标量的湍流、波动的影响，可以使用不同类型的封闭模型。与单相流动相比，多相流动动量方程中所模拟的项数是非常大的，这使得多相流模拟中的湍流模型非常复杂。

在k－ε模型内常用三种方法模拟多相流中的湍流：

（1）混合湍流模型；

（2）分散湍流模型；

（3）每相湍流模型。

模型的选择依赖于应用中第二相湍流的重要性。

注意：下面给出的每一种方法的描述都是基于标准k－ε模型的。将多相流的湍流修正为基于重整化群（RNG）k－ε和可实现（realizable）k－ε模型时可作相似的处理，因此这里不再明确地给出。

1）混合湍流模型

混合湍流模型是默认的多相湍流模型。它代表了单相k－ε模型的第一扩展，应用于相分离、分层（或接近分层）的多相流以及相之间的密度比接近于1的情况。这种情形下，使用混合属性和混合速度捕获湍流的重要特征是足够的。

描述这个模型的k－ε方程如下：

和

这里，混合密度ρm和混合速度v　→m由下式计算：

湍流粘度μt，m由下式计算：

湍流动能的产生项Gk，m由下式计算：

这些方程中的常数与单相k－ε模型中的描述相同。

2）分散湍流模型

当第二相为稀相时，分散湍流模型是合适的模型。这种情形下，颗粒间的碰撞可忽略，而对第二相随机运动起支配作用的是主相湍流的影响。所以第二相的波动量根据主相的平均特征、颗粒松弛时间和粒子相互作用时间给出。

模拟湍流的分散方法涉及以下假设：

①采用修正的k－ε模型描述连续相：连续相湍流预测是通过使用标准k－ε模型并补充包含相间湍流动量传递的附加项获得的。

②采用Tchen－theory关系式描述离散相：离散相湍流度的预测是采用Tchen的均匀湍流离散粒子分量理论来获得的。

③相间湍流动量传递：在湍流多相流动中，动量交换项包含了离散相瞬态分布和湍流流体运动之间的关系。需要考虑湍流流动所引起的离散相输送的分散问题。

④相加权平均方法：在模拟湍流多相流的分散时平均方法的选择是有影响的，两步平均法会导致出现相体积分数的波动。然而，当使用两步平均法加上对湍流的相加权平均时，体积分数的湍流波动不会出现。

（1）连续相中的湍流

涡粘性模型常用于计算平均波动量。连续相q的雷诺应力张量采用如下形式：

这里，是相加权速度。

根据q相的湍流动能，湍流粘度μt，q写成如下形式：

载能湍流涡的特征时间定义如下：

这里，εq是耗散率，Cμ＝0.09。

湍流涡的长度标尺为

湍流预测从修正的k－ε模型获得

和

这里，代表了分散相对连续相q的影响

代表了湍流动能的产生，所有其他项与单相k－ε模型中的项有相同的意义。

项可从连续相的瞬态方程获得并采用如下形式，其中M代表第二相的数量

它可以简化为

这里，kpq是连续相q的分散相p的速度的协方差［由后面的方程（6－76）计算而得］是相对速度是漂移速度。

根据Elgobashi等人提出的模型：

这里，C3ε＝1.2。

（2）离散相中的湍流

表征运动的时间和长度标尺用于估计分散（dispersion）系数，相关函数和每一分散相的湍流动能。和作用于离散相p的惯性影响相连接的特征粒子的松弛时间定义为

沿着颗粒轨道计算的拉格朗日（Lagrangian）积分时间标尺主要受交叉轨道的影响，其定义为

式中

Cβ＝1.8－1.35cos2θ　　（6－73）

其中，θ是平均颗粒速度和平均相对速度的夹角。这些方程中的两个特征时间比值写做

根据Simonin的研究，可将离散相p的湍流量度写为如下形式：

式中，CV＝0.5是附加的质量系数。

（3）相间湍流动量传递

多相流的湍流曳力项［方程（6－11）中的Kpq（v　→p－v　→q）］表达式如下，对于离散相p和连续相q：

方程（6－80）右边的第二项包含漂移速度

式中，Dp和Dq是扩散率，σpq是湍流施密特（Schmidt）数。当在多相流中使用Tchen的理论时，可假设Dp＝Dq＝Dt，pq，并默认σpq的值为0.67。

在体积分数中漂移速度起因于湍流波动，当乘以交换系数Kpq时，它用做湍流动量交换项的修正。

3）每相湍流模型

最普通的多相湍流模型为每一相求解一套k－ε输运方程。当湍流传递在相间起重要作用时，这个湍流模型是合适的选择。

（1）输运方程

雷诺应力张量和湍流粘度采用方程（6－61）和（6－62）计算。湍流预测从下式获得：

Clq和Cql可以近似地写为

这里，ηlq由方程（6－74）定义（只是用l代替p）。

（2）相间湍流动量传递

湍流曳力项［方程（6－11）中的如下，其中l是分散相［代替方程（6－11）中的p］，q是连续相：

这里是相加权速度，v　→dr，lq是相l的漂移速度［使用方程（6－81）计算，用l取代p］。

如前文指出的，漂移速度起因于体积分数的湍流波动。当乘以交换系数Klq时，用于修正湍流中的动量交换项。

6.3　VOF模型

6.3.1　模型的应用背景及特点

所谓VOF模型，是一种在固定的欧拉网格下的表面跟踪方法，该模型的构建基于欧拉框架。当需要得到一种或多种互不相融流体间的交界面时，可以采用这种模型。在VOF模型中，不同的流体组分共用着一套动量方程，计算时在全流场的每个计算单元内都记录下各流体组分所占有的体积率。

VOF模型的应用例子包括分层流、自由面流动、灌注、晃动、液体中大气泡的流动、水坝决堤时的水流、对喷射衰竭（jet breakup）（表面张力）的预测，以及求得任意液－气分界面的稳态或瞬时分界面。该模型能够比较好地反映多相流之间的界面情况，比如大的气泡以比较慢的速度在液体中流动、气－液界面等。由于VOF模型采用的方程中的各项物性参数，如密度、粘度等，是各相物性的体积平均值，所以要求各相的速度之间差别不能太大，否则会对计算结果精度的影响很大。一般情况下VOF采用非稳态模拟比较好。主相的体积值不是从体积守恒方程得到的，而是1减去其他离散相的值。

6.3.2　主要数学模型

VOF模型依赖的是两种或多种流体（或相）没有互相穿插（interpenetrating）这一事实。对拟增加到模型里的每一附加相引进一个变量，即计算单元里的相的体积分数（volume fraction of the phase）。在每个控制容积内，所有相的体积分数的和为1。所有变量及其属性的区域被各相共享并且代表了体积平均值（volume－averaged　values），只要每一相的体积分数在每一位置是可知的。这样，在任何给定单元内的变量及其属性或者纯粹代表了一相，或者代表了相的混合，这取决于体积分数值。换句话说，在单元中，如果第q相流体的体积分数记为αq，那么下面的三个条件是可能的：

①αq＝0：第q相流体在单元中是空的。

②αq＝1：第q相流体在单元中是充满的。

③0＜αq＜1：单元中包含了第q相流体和一相或者其他多相流体的界面。

基于αq的局部值，适当的属性和变量在一定范围内分配给每一控制容积。

6.3.2.1　守恒方程

1）连续性方程

在VOF模型中，跟踪相之间的界面是通过求解一相或多相的体积分数的连续性方程来完成的。对第q相，这个方程如下：

式中是第p相到第q相的质量输送则是第q相到第p相的质量输送。一般情况下，方程（6－86）右端的源项为零。体积分数方程不是为主相求解的，主相体积分数的计算基于如下的约束：

在两相流系统中，如果相用下标1和2表示，第二相的体积分数会被跟踪，那么每一单元中的密度由下式给出：

通常，对n相系统，体积分数平均密度采用如下形式：

所有的其他属性［例如粘度（viscosity）］都以这种方式计算。

2）动量守恒方程

通过求解整个区域内的单一动量方程，作为结果的速度场是由各相共享的。如式（6－90）所示，动量方程取决于通过属性ρ和μ所表达的各相的体积分数

近似共享区域的一个局限是，当各相之间存在大的速度差异时，靠近界面处的速度的计算精确性会受到不利影响。

3）能量守恒方程

能量方程，也是各相共享的，表示如下

VOF模型处理能量E和温度T的方法与处理质量平均变量一样：

式中每一相的Eq是基于该相的比热和共享温度。属性ρ和keff（有效热传导）是被各相共享的；源项Sh包含辐射的贡献，也有其他体积热源。和速度场一样，在相间存在大的温度差时，靠近界面的温度的精确度也受到限制。在属性有几个数量级的变化时，这样的问题也会增长。例如，如果一个模型包括液体金属和空气，材料的导热性有4个数量级的差异，如此大的差异使得方程需要设置各向异性的系数，以达到收敛性和精度范围。

6.3.2.2　界面附近的插值

控制体积公式要求计算穿过控制体积面的对流和扩散通量，并与控制体积本身内部的源项平衡。对VOF模型，有4种常用的方案计算面的通量：几何重建（geometric reconstruction）方案、物质接受（donor－acceptor）方案、欧拉显式方案和欧拉隐式方案。

在几何重建和物质接受方案中，用了特殊的插值处理两相之间界面附近的单元。图6－1显示了用这两种方案计算的过程中，沿着假定的界面的实际界面的形状。

图6－1　几何重建和物质接受方案的界面形状

欧拉显式和隐式方案以相同的插值方式处理这些完全充满一相或其他相的单元（也就是，使用标准迎风、二阶迎风或者QUICK方案），而不采用特殊的处理。

1）几何重建方案

在几何重建方案中所使用的标准插值方法用于获得界面通量，无论何时，一旦单元被某一相充满，标准插值方案将被用于获取界面通量。当单元靠近两相之间的界面时，应使用几何重建方案。

几何重建方案使用分段线性的方法描绘了流体之间的界面。这个方案是最精确的并适用于通用的非结构化网格。几何重建方案是从Youngs关于非结构化网格的研究工作中归纳出来的。它假定两流体之间的界面在每个单元内有个线性斜面，并使用这个线性形状来计算穿过单元面的流体的对流。

几何重建方案的第一步是根据体积分数和单元引出的信息计算相对于每个部分充满单元的中心的线性界面的位置；第二步是使用已计算的关于面上的法向和切向速度分布的信息线性的界面描述来计算穿过每个面的流体的对流量；第三步是使用前面的步骤中计算的通量平衡计算每个单元的体积分数。

当使用几何重建方案时，必须计算时间依赖。同样，如果使用正投影网格（也就是如果网格节点的位置在两个子区域相交的边界上是一样的），必须确保在区域内没有双边（零厚度）的壁面。如果有，必须分开它们。

2）物质接受方案

在物质接受方案中，使用的标准插值方法用于获得面的通量，无论何时单元内完全充满一相或其他相。当单元靠近两相之间的界面时，物质接受方案用于决定穿过面的流体的水平对流量。这个方案把一个单元看作一定数量的流体来自一相的捐赠（donor），把相邻的单元看作相同数量流体的接受（acceptor），这样使用防止了界面上的数值扩散。来自对流跨过一个单元边界一相流体的数量受限于两个值的最小值：捐赠单元的充满体积和接受单元的自由体积。

界面的方向也用于决定面的通量。界面的方向是水平的还是垂直的，取决于单元内第q相的体积分数梯度的方向和问题中共享面的相邻单元。依靠界面的方向和它的运动，通过纯的迎风、纯的顺风或二者的联合获得通量值。

当使用物质接受方案时，必须计算时间依赖。而且，物质接受方案仅用于四边形和六面体网格。另外，如果使用了正投影网格（也就是如果网格节点的位置是一样的在两个子区域相交的边界上），必须确保在区域内没有双边（零厚度）的壁面。如果有，必须分开它们。

3）欧拉显式方案

在欧拉显式方案中，标准的有限差分插值方法被用于前一时间步的体积分数的计算

式中，n＋1为新时间步的指标；n为前一时间步的指标；αq，f为一阶迎风、二阶迎风、QUICK算法中的第q相体积分数的面值；V为单元的体积；Uf为以法向速度通过面的体积流量。

这个公式在每一时间步上不需要输送方程的迭代解，在隐式方案中是需要的。同样，当使用欧拉显式方案时，必须计算时间依赖。

4）欧拉隐式方案

在欧拉隐式方案中，标准的有限差分插值方法用于获得所有单元的面通量，包括那些界面附近的

由于这个方程需要当前时间步的体积分数值（而不是上一时间步，如欧拉显式方案），在每一时间步内标准的标量输送方程为每一个第二相的体积分数迭代性地求解。

欧拉隐式方案可用于时间依赖和稳态的计算。

6.3.2.3　表面张力和壁面粘附

VOF模型也可以模拟沿着每一对相之间的表面张力的影响，主要是通过附加的描述相和壁面之间的接触角来实现。

1）表面张力

由于流体中分子之间的引力，产生了表面张力。例如，考虑水中的一个气泡。在气泡内，由于其周围相邻分子的作用，作用在分子上的净力为零。然而，在表面上，净力是呈放射状地向内的，跨过整个球面的径向分力的联合影响是表面收缩，因而增强了表面凹侧的压力。表面张力是一种仅作用在表面上的力，在这个例子中它必须是保持平衡的，它扮演了平衡内部放射状的分子引力和跨过表面的放射状的外部压力梯度的角色。在两种流体分离的地区（但是它们当中的一个不是球泡的形式），表面张力的作用是通过减小界面的面积使自由能最小化。

常用表面张力模型是由Brackbill等人提出的连续表面张力模型。用这个模型，VOF计算中附加的表面张力导致了动量方程中的源项。为了理解这个源项的起源，考虑沿着表面张力为常数的特殊情况，那些地方只考虑垂直于界面的力。可以看出，跨过表面的压降依赖于表面张力系数σ以及通过两个半径的正交方向量度的表面曲率R1和R2

这里p1和p2是两种流体界面两侧的压力。

使用连续表面张力（CSF）模型公式时，这里的表面曲率是从垂直于界面表面的局部梯度计算的。n为表面法线，定义αq为第q相体积分数的梯度

n＝▽αq　　　　（6－96）

表面曲率κ是为了区别单位法向量＾n而定义的

κ＝▽·＾n　　　　（6－97）

式中

表面张力也可以根据越过表面的压力跳跃写出。表面张力使用散度定理可以表示为体积力，正是这个体积力成了添加给动量方程的源项。它有如下形式：

该式允许在多于两相存在的单元附近光滑地叠加。如果一个单元中只有两相，那么κi＝－κj，▽αi＝－▽αj，式（6－99）可简化为

这里，ρ是使用式（6－99）计算出的体积平均密度。式（6－100）显示出一个单元表面张力源项是与单元的平均密度成比例的。

注意，三角形和四面体网格上表面张力影响的计算不如四边形和六面体网格的计算精确，所以表面张力影响最重要的地区应当采用四边形和六面体网格。

表面张力影响重要性的决定因素是两个无量纲数：雷诺数Re和毛细数（capillary number）Ca或雷诺数Re和韦伯数（Weber number）We。当Re《1时，起主要作用的是毛细数

当Re》1时，起主要作用的是韦伯数

这里，U是自由流速度。如果Ca》1或者We《1，表面张力效应可以忽略。

2）壁面粘附

与表面张力模型联合时选择指定一个壁面粘附角在VOF模型中也是有用的。这个模型假定流体与壁面产生的接触角常用于调整壁面附近单元表面的法向，而不是加强壁面本身的边界条件。这个所谓的动力壁面边界条件导致了壁面附近表面曲率的调整。

如果θw是壁面的接触角，那么挨着壁面的实际单元的表面法向为

这里，＾nw和t　＾w分别是壁面的单位法向量和切向量。这个接触角与非壁面处单元通常计算表面法向的联合决定了表面的局部曲率，这个曲率常用于调整表面张力计算中的体积力项。

接触角θw为壁面和壁面上界面切线的夹角，如图6－2所示。

图6－2　接触角示意图

6.4　颗粒轨道模型

颗粒轨道模型基于欧拉－拉格朗日框架，该模型认为流体相是连续的，颗粒相是离散的。因此，相对于欧拉模型中固相连续性的假设，欧拉－拉格朗日模型对颗粒流体系统的模拟有更合理的理论基础。

根据气固两相系统中颗粒运动的特点，该模型将两相流动中颗粒的运动过程分解为受冲力支配的瞬时碰撞运动及受流体曳力控制的悬浮运动，从而建立了颗粒运动分解模型。该模型中，流体的运动规律用连续介质的N－S方程进行描述，而颗粒的行为则通过在拉格朗日坐标系中分析每一个单颗粒的运动轨道而进行描述。其中，在颗粒与颗粒相互作用的过程中，其运动规律服从碰撞动力学中的动量守恒定律；在流体与颗粒相互作用的悬浮过程中，颗粒在曳力、重力等力的作用下，其运动规律服从牛顿动力学中的力平衡方程。这样，每个颗粒的速度及位移的更新由邻近颗粒对它的碰撞过程及流体对它的悬浮过程来确定。与此同时，颗粒对流体的瞬时作用则反映在离散的、两相耦合的N－S方程不断得以修正的数值求解过程中。

该模型的基本方程请参阅本书5.4节的内容，下面将针对颗粒轨道模型中几种常用的处理颗粒与颗粒之间相互作用的数学模型进行详细介绍。

6.5　描述颗粒之间相互作用的模型

当采用拉格朗日方法对固相运动进行描述时，随着固相浓度的增加，颗粒间的碰撞对气固两相流动有非常重要的影响，因此必须考虑颗粒与颗粒之间的作用。

根据处理颗粒碰撞方式和计算方法的不同，现阶段计算机模拟主要采用两种模型，即硬球模型和软球模型，将在下文进行详细介绍。

此外，对于如何判断颗粒是否发生碰撞的问题，主要有两种解决方法，即确定性方法和随机性方法。

确定性方法的基本思想是：当两个运动颗粒间的距离小于一定值时将发生碰撞，碰撞后颗粒运动速度按经典力学的规律计算得到。这种方法的优点在于只要颗粒的初始状态给定，就能够计算出所有颗粒的运动。但是采用这种方法模拟任意一个颗粒的运动时，都要考虑到所有其他颗粒的运动以发现是否有可能与其发生碰撞，于是模拟所需的计算时间很长，因此通常确定性方法所用的颗粒数都不能太大，使得确定性方法的应用范围受到了限制。

随机性方法的基本思想是：将颗粒间的相互碰撞和气体分子运动论中的分子间的碰撞相类比，认为颗粒的碰撞具有随机性。用于颗粒碰撞处理的随机性方法目前一般采用直接模拟蒙特卡洛（Direct　Simulation　Monte　Carlo，DSMC）方法。该算法用相对少量的取样颗粒代替大量的真实颗粒进行模拟，使得模拟计算效率大为提高。

下面将针对在拉格朗日框架下的三种描述颗粒间相互作用的模型或方法，即软球模型、硬球模型和DSMC方法进行详细介绍。

6.5.1　软球模型

6.5.1.1　模型的应用背景及特点

软球模型是目前小规模数值计算中最接近实际的一种模型。它从最基本的物理定律出发来考虑颗粒间的受力，并通过运动轨迹跟踪流场内的每一个颗粒，能够比较准确、细腻地描述颗粒碰撞接触过程的作用力，其模拟计算的固相浓度能达到空间饱和（堆积）状态，因而适用于密相流化床、密相气力输送和料仓中重力流动等稠密气固两相流动的建模与计算。该模型一个很大的优点是可以处理多个颗粒同时碰撞，这在模拟准静态系统时（如固定床）是非常重要的。但该模型的计算量大，所模拟颗粒数量、颗粒大小等参数受到严重限制，特别是针对涉及超浓相、粉体颗粒的数值模拟时，计算量更是巨大，难以进行接近实际规模的数值模拟计算。因此，目前软球模型的应用成果仍局限于大尺寸颗粒、小规模装置以及短时间的模拟，尚不能用于工程意义上的稠密气固两相流问题。

软球模型最初是根据Cundall等提出的应用于土壤力学的离散单元法（Distinct Element Method，DEM）发展而来的，其基本思想是：假设颗粒碰撞可以持续一定的时间并允许颗粒在碰撞的过程中发生轻微的重叠现象，采用弹簧、阻尼器和活塞来模拟颗粒与颗粒之间和颗粒与壁面之间的碰撞过程，如图6－3所示。

6.5.1.2　主要数学模型

软球模型通过跟踪离散颗粒场中的每一个颗粒来判断颗粒是否发生碰撞，当两个球形颗粒发生弹性对心碰撞时，首先在接触点发生弹性变形，颗粒在前进方向上受到阻力，该阻力的大小与法向变形位移和材料的刚度成正比，在达到最大变形位移时，颗粒停止运动，随后在该力作用下，沿原来的运动直线反弹。

图6－3　模拟接触力的力学系统

对于非完全弹性碰撞，碰撞后颗粒的动能发生损失，动能损失的大小与颗粒材料的物性和碰撞时的相对速度有关，该部分动能损失可归纳为在碰撞过程中受到一个与颗粒方向相反的反力，该力的大小等于两颗粒相对速度与阻尼系数的乘积，当两颗粒发生偏心碰撞时，相撞点处的接触力可分为法向分力和切向分力，分别由法向变形位移和切向变形位移以及法向动能损失和切向动能损失进行计算。法向分力的作用结果如同对心碰撞，切向分力的作用结果是对颗粒球心产生一个矩，该矩使颗粒发生旋转，由该矩和颗粒的转动惯量可求出所产生的角加速度。

切向力的极值受到颗粒表面摩擦系数与法向分力乘积的限制，当所计算出的切向分力大于该乘积时，两颗粒在接触表面将发生滑移。更一般的情况是两个旋转颗粒发生偏心碰撞。这时除计算法向位移和切向位移外，还应计及由于颗粒自转在接触点处所造成的切向速度。当一个颗粒同时与几个颗粒相撞时，通过矢量叠加可以算出颗粒所受到的合力与合力矩。

上述物理过程可通过以下数学模型描述：

当颗粒i与颗粒j的质心间距小于两颗粒半径之和时，则发生碰撞。在碰撞点处产生力的作用，记为颗粒间碰撞力fCij，该力可以分为法向力fCnij和切向力fCtij，法向分力使颗粒发生平动，而切向分力产生矩的作用，使颗粒发生转动

式中，δnij和δtij分别为颗粒法向和切向的变形量；vrij为颗粒的相对速度；vsij为颗粒碰撞点的滑移速度；kn和kt分别为法向和切向的弹性系数；ηn和ηt分别为法向和切向的阻尼系数，可由恢复系数e等确定；nij为单位法向量；ωi和ωj分别为两颗粒的旋转角速度。当｜fCtij｜＞μ｜fCnij｜时，颗粒相对滑移，滑移产生的切向分力为

式中，tij＝vsij/｜vsij｜，μ为滑动摩擦系数。

F＝fC＋fF　　　　（6－109）

v＝F/m＋g　　　　（6－110）

ω＝T/I　　　　　（6－111）

式中，F为合力，fF为流体力，m为颗粒质量，g为重力加速度，T为合力矩，I为颗粒转动惯量。

v＝v0＋6vΔt　　　　（6－112）

r＝r0＋6rΔt　　　　（6－113）

ω＝ω0＋6ωΔt　　　　（6－114）

式中，v为颗粒速度向量，Δt为时间步长，r为颗粒重心的位置向量，下标0为前一Δt旧值。

更为常见的是颗粒i同时与n个颗粒发生碰撞，这时只要对多个颗粒中的每一个颗粒j分别计算与颗粒i的碰撞受力，然后计算作用于颗粒i的合力和合力矩。

6.5.1.3　模型参数的确定

根据颗粒间接触作用力模型，包含了摩擦系数、弹性系数、阻尼系数、时间步长等参数，其参数的确定如下：

（1）摩擦系数

颗粒表面之间的摩擦系数f与材料的物理性质及颗粒表面的特性有关。摩擦系数的求解公式为

f＝tanφ　　　　（6－115）

其中，φ为颗粒休止角。

（2）弹性系数

颗粒的弹性系统是由弹性模量（E）与泊松（Possion）比（v）决定的。根据Hertz的理论，得出其法向弹性系统的求解公式为

其中，

根据Mindlin等人的理论，其切向的弹性系数为

其中，δij为两颗粒间的法向位移，Gs为剪切模量

表6－1中列出了各文献中k值的选取情况。

表6－1　文献中k值的选取

（3）阻尼系数

如果一个弹性系统中有了动能，则系统会在平衡位置附近作间谐振动。为了避免这一点，就要采用加阻尼的办法来耗散系统在振动过程中的动能。特别是在静力问题中，必须加阻尼来吸收系统中的动能，以使系统达到稳定的状态。颗粒与颗粒及颗粒与壁面间碰撞时的阻尼系数由表6－2给出。其中，ηn，ηs分别表示法向、切向阻尼系数，e为碰撞弹性恢复系数，定义为

其中，v、v′分别为颗粒碰撞前后的速度。

表6－2　颗粒之间及颗粒与壁面之间碰撞的阻尼系数

通常颗粒的切向碰撞恢复系数与法向碰撞恢复系数相等。对于理想情况，碰撞恢复系数为1，即没有机械能损失。对于给定颗粒物性的情况，碰撞恢复系数是一个与颗粒质量、弹性系数和阻尼系数无关的常数值。

需要注意的是，颗粒离散单元法中引入阻尼系数是为了提供耗能装置，并非为了得到准静态解，因此阻尼系数的选取应具有一定的灵活性，以满足最大程度模拟为原则。

（4）时间步长

在软球方法中，时间步长Δt是影响数值计算的一个关键因素，Δt太小，会使得计算量很大；Δt太大，会使得计算发散。

Tsuji等人根据式（6－122）确定时间步长，它与颗粒刚度系数k和颗粒质量m有关

K.D.Kaufi等人根据式（6－122）确定时间步长

其中，λ＝0.163　1v＋0.876　6，v为泊松比，G为剪切（shear）模量，ρp为颗粒密度，dpmin为颗粒的最小直径。

6.5.2　硬球模型

6.5.2.1　模型的应用背景及特点

硬球模型基于冲量守恒定律，认为颗粒间碰撞是瞬时发生的。与软球模型相比，硬球模型不考虑颗粒受力变形等细节，而是采用恢复系数和摩擦系数直接求解颗粒碰撞后的速度，并假设颗粒的碰撞过程是在一个时间步长内完成的，因而有效地节省了计算时间以及计算机内存。然而，该模型假设颗粒间的碰撞为两两碰撞，所以无法处理三体或三体以上同时碰撞的情况。因此，硬球模型仅适用于低浓度、快速颗粒流的模拟，而不能针对稠密气固两相流动进行计算。

6.5.2.2　主要数学模型

设两个颗粒a，b，质量分别为ma，mb，半径分别为ra，rb，碰撞前的线速度和角速度分别为碰撞后的线速度和角速度分别为va，ωa，vb，ωb，而e，ewall分别为颗粒与颗粒、颗粒与壁面碰撞时的恢复系数。

在两个颗粒碰撞时做以下假设：

（1）颗粒是球形的，且是准刚性的，碰撞后的颗粒不变形。（2）假设颗粒碰撞是二体的瞬时碰撞。

（3）碰撞仅发生在两个颗粒的接触点上。

（4）颗粒间的相互作用是瞬时冲力，在碰撞过程中忽略其他力。

硬球模型中，颗粒间的碰撞是二体瞬时碰撞，颗粒间的碰撞满足动量定理

其中，n为碰撞时颗粒a的质心指向颗粒b的质心的单位法向量；J是施加在颗粒上的冲量；I＝0.4mr2，是转动惯量和碰撞前颗粒的位置是给定的，只有冲量J和碰撞后的颗粒速度va，ωa，vb，ωb未知的。

为了求解上述方程的未知量，再做以下假设：

（1）在颗粒碰撞时，颗粒质心间的距离为两颗粒的半径和，且颗粒不发生变形。

（2）颗粒的滑动摩擦服从库仑定律。

（3）一旦一个颗粒停止了滑移，就不会有进一步的滑移。

设和G分别为颗粒碰撞前和碰撞后的相对速度，f为摩擦系数，t为颗粒相对速度的切向单位向量是颗粒碰撞前相对速度的切向分量。

当两个颗粒在碰撞时有滑移时，碰撞后的速度公式为

当两个颗粒在碰撞时停止滑移，碰撞后的速度公式为

其中，Gct＝G－（G·n）n＋riωi×n＋rjωj×n。

判断颗粒碰撞时是否滑移的条件是

判断颗粒碰撞时停止滑移的条件是

6.5.3　直接模拟蒙特卡洛（DSMC）方法

6.5.3.1　模型的应用背景及特点

DSMC方法是在硬球模型的基础上结合了判断颗粒碰撞的蒙特卡洛方法。传统方法是根据颗粒所处位置是否有重叠的部分来判断颗粒之间是否发生碰撞，而DSMC方法以样本颗粒代替真实颗粒，并通过概率来判断碰撞。该方法的计算工作量较软球模型、硬球模型大为减少，因而有效地节省了计算时间与计算机内存，可对大规模颗粒构成的稀疏气固两相流动进行数值模拟。但由于该方法基于硬球模型，因而不能针对稠密气固两相流动进行计算。

6.5.3.2　主要数学模型

1）DSMC方法的思想

DSMC方法的基本要点可以简述如下：

用少量的样本颗粒代替众多真实颗粒，通过统计样本颗粒的运动状态实现对真实颗粒运动的模拟。这里，为了用少数样本颗粒来代替众多真实颗粒重现碰撞，如图6－4所示，用一个简单的样本颗粒场来代替真实颗粒场。图6－4中，（a）是在一个网格中的真实颗粒场，这个场中每个颗粒都有不同的速度，令颗粒数目为n；（b）表明在同一网格中具有不同速度的三个样本颗粒，分别用白、灰、黑三种颜色描述，令样本颗粒数目为N；在样本颗粒场（c）中，同一网格中每种样本颗粒（具有相同速度）是随机分布的，每一种样本颗粒所代表的真实颗粒数目为n1＝n/N，n2＝n/N，…DSMC方法的基本思想就是用颗粒场（c）来代替颗粒场（a）。可以用数学方法证明如果样本颗粒的数目充足大时（c）与（a）在统计意义下相同。一旦真实场用样本场代替，碰撞概率理论就可以用于颗粒间的碰撞。因此，在颗粒场（b）的基础上来计算是可行的。在DSMC方法中，轨道计算可应用于样本颗粒的计算。对于颗粒场（a）的情况，只要计算（b）中的三个颗粒即可。一个样本颗粒是否与另一个碰撞，完全由颗粒的碰撞概率决定。

图6－4　真实颗粒场与样本颗粒场示意图

2）颗粒的碰撞概率

在DSMC方法中，判断一个颗粒是否与另一个颗粒碰撞主要由碰撞概率来决定，而不是像确定性方法那样由颗粒的轨道来决定。

现有处理颗粒碰撞的概率模型主要有两类：一类是基于气固两相流中颗粒间的碰撞，如Kitron模型；另一类模型是在拉格朗日轨道的基础上建立的颗粒间的相互碰撞模型，如Tanaka模型、Oesterle模型和Sommerfeld模型。

其中第一类模型的相似假设限制了其适用范围，一般仅适用于高雷诺数、稀相、大颗粒（直径＞100μm）的气固两相流动系统，而第二类模型既可适用于稀相，也可用于浓相，而且不受颗粒以及雷诺数大小的限制。

在拉格朗日轨道的基础上建立的颗粒随机碰撞模型的具体描述如下：

（1）Sommerfeld等提出的模型：在体积为V，颗粒数为N的网格中，在t时间内任何一对颗粒发生碰撞的频率为

其中，ri，rj为相碰颗粒的半径；｜Gij｜为这一对颗粒相对速度的模。样本颗粒的速度服从正态分布，其方差和均值采用统计足够多的、在一定时间内经过所给网格的颗粒速度的方法求得。

（2）Oesterle等提出的模型：半径为ri的颗粒i与环境中其他颗粒发生碰撞的频率为

其中，系数是把一维情况扩展到三维情况的修正系数。假设颗粒的相对速度关于局部平均颗粒速度服从Maxwellian分布，样本颗粒的速度及旋转速度取该网格的平均速度，两颗粒的碰撞点采用随机方法求得。

（3）Tanaka提出的模型

其中，nσ为颗粒浓度，Gij为其相对速度，Dp是颗粒的直径。

因为Sommerfeld模型和Oesterle模型必须满足所要求的分布，因此操作起来比较复杂，而Tanaka模型相对简单易行。Tanaka模型具体的物理意义可描述如下：

假设系统中有n个真实颗粒，用N个样本颗粒来代替，且颗粒被假设为大小均匀、刚性且完全光滑的球体，颗粒的直径为Dp。设颗粒i具有速度vi，而一组颗粒j都有相同的速度vj，且在空间随机分布。这种情形与颗粒i有速度Gij＝vi－vj，而颗粒j静止的情形相同。颗粒i在Δt内移动的距离是GijΔt。

这样，颗粒i与颗粒j碰撞的概率即可由图6－5表示，其表达式为

图6－5　Tanaka模型的颗粒碰撞概率示意图

其中，nσ是一个样本颗粒代表真实颗粒数目的比值，即

；V代表流场网格的体积，即V＝Δx×Δy×Dp。故Δx和Δy分别代表将流场分隔成小网格的宽度和高度，由于是二维运算，所以在这里以一个颗粒的直径为虚拟厚度。

3）DSMC方法的基本假设

DSMC方法的基本假设如下：

（1）将流场区域划分成小的网格区域，碰撞概率的计算和碰撞对的寻找都在这些小网格内进行。

（2）对实际颗粒进行抽样，计算的样本颗粒数目远小于实际颗粒数目。

（3）颗粒的碰撞时间间隔Δt远小于颗粒的悬浮运动时间间隔，以使得颗粒间的相互碰撞和颗粒的自由悬浮运动解耦。

4）修正Nanbu算法

（1）修正Nanbu算法的原理

用于搜索颗粒的碰撞对的蒙特卡洛方法有许多种，如时间统计法、Nanbu算法和修正Nanbu算法。和别的方法相比，修正Nanbu算法只产生一次均匀随机数，这样就节约了系统资源，简化了计算步骤，因此选用修正Nanbu算法来搜索颗粒碰撞对。

图6－6为颗粒碰撞概率示意图，图中将一个单元长度的线段分成N等份，即每个线段的长度为1/N，其中N是典型颗粒的数目。每一小线段的长度用来与颗粒i和每个典型颗粒碰撞的概率相比较。由式（6－141）知，颗粒i与颗粒1，2，…，j，…与N的概率为Pi1，Pi2，…，Pij，…，和PiN，分别表示为

图6－6　颗粒碰撞概率示意图

颗粒i与其他颗粒的碰撞概率为以上概率之和，即为

这个概率Pi满足

0＜Pi＜1

如果时间步长Δt足够小，以上概率Pi1，Pi2，…，Pij，…，和PiN总满足条件

根据以上描述，取R为在［0，1］均匀分布的随机数，如果，我们就只计算Pij；如果有，则认为颗粒i与颗粒j不碰撞；反之，如果，则认为颗粒i与颗粒j发生碰撞。这里要求所有的可以通过选择足够小的Δt来满足这个要求，这个算法的计算时间为O（N）。

（2）修正Nanbu算法的描述

详细的算法可以这样描述：产生一个［0，1］的均匀分布随机数列，从中取出一个随机数R，一个候选碰撞颗粒的编号由下式给出

j＝int［R×N］＋1　　　　（6－144）

若

则满足颗粒i与颗粒j在Δt内发生碰撞的要求。

5）随机数发生器

随机数是蒙特卡洛方法的基础。模拟中需要产生在［0，1］之间的均匀随机数。杨自强等人提供了各种随机数发生器，它们都能产生［0，1］之间的均匀随机数，其中简单可行的有：

（1）线性同余法随机数发生器，其递推式为

线性同余法历史悠久，算法简单，给定模M后，其余参数的选择应保证｛Xi｝有最大的周期。但是采用线性同余法时有长周期相关现象。

（2）进位加随机数发生器，其递推式为

进位加随机数发生器的随机数产生速度快，并且有惊人的长周期，这是因为递推式里没有乘法运算。但这种发生器可以理解为模M很大时的线性同余法随机数发生器。

（3）混沌映射随机数发生器，其递推式为

混沌映射随机数发生器的序列｛yi｝是无限不循环的，所以它可以产生周期为无限的U（0，1）分布的随机数列。这个方法等价于a＝1且c＝1－Ln2/Ln3＝0.369　070…的线性同余法随机数发生器，但实际应用中，c只能取有限位，从而得到的序列不可能无限不循环，而因子a＝1，显然这种发生器的统计结果不可能优于两个参数（a，c）都经过选择的线性同余法随机数发生器。

（4）线性同余组合随机数发生器，其递推式为

线性同余组合随机数发生器优于组成它的任一随机数发生器，但是它本质上等价或近似等价为一个单个的关联线性同余随机数发生器，不过它可以具有巨大的模其中m表示线性同余随机数发生器的个数，用M（j）表示第j个随机数发生器的模。