CN109840375B - 一种液固流化床cfd曳力模型的确认方法 - Google Patents

Abstract

本发明公开了涉及流化床多相流CFD模拟领域的一种液固流化床CFD曳力模型的确认方法。基于对颗粒群体系中单颗粒的受力分析建立力学方程，结合CFD曳力模型表达式，求出床内颗粒体积浓度，并将其与实验测量值进行比较，最终确认出适用于该物性参数和操作条件下的最佳CFD曳力模型。该方法通过理论计算替代CFD计算，在确保准确性的前提下，克服了当前液固流化床CFD模拟中曳力模型的确认耗费大量时间与计算资源的弊端；同时也实现了对多个CFD曳力模型的确认。

Description

一种液固流化床CFD曳力模型的确认方法

技术领域

本发明涉及流化床多相流CFD(Computational Fluid Dynamics)模拟领域，公开了一种液固流化床CFD曳力模型的确认方法。

背景技术

由于具有良好的颗粒混合特性以及较高的传热传质效率，液固流化床的应用几乎遍及整个过程工业，包括了化工、能源、矿业、环保以及生化等领域。在这些应用中，颗粒流体多相流显著影响着床中的传热传质与化学反应等过程，其居于基础地位；因此，深刻理解床内流体动力学特性对液固流化床的操作、设计、优化以及放大起到至关重要的作用。

液固床内流体动力学的研究主要分为实验和CFD模拟两类。实验法是对物理真实的直观记录，获取数据可信度高，但其成本巨大，耗费大量的时间与人力、物力，而且实验测量存在较大的时空局限，往往难以获得全面而细致的流场信息。

CFD技术是通过有限差分或者有限体积法将描述流体运动的偏微分控制方程离散后，利用计算机进行数值求解，得到压力、速度等信息，实现流场的再现，从而揭示出流体的流动特征和演化规律。与实验研究相比较，CFD模拟的固有优势在于其不受实验测量技术的限制，可以高效而低成本地获得宏观流动特征与微观流场细节。随着计算机技术的发展与数值算法和物理模型的日益完善，CFD已经从作为实验方法的有益补充，逐渐发展为流体动力学研究的重要手段和工程实践中不可或缺的工具。

尽管CFD技术业已受到广泛认可，但其本质上仍然是对物理真实的虚拟再现，所以计算结果的可信度就成为关乎CFD模拟成败的头等大事。然而令人遗憾的是，CFD使用者对计算结果的可信度并不能轻易给出答案，只有经过专业的严格的验证和确认后，计算结果方可得到认可，亦即CFD模拟是可信的。

根据美国航空航天学会CFD标准委员会的定义，可信度研究的基本内容和方法就是CFD的验证和确认。CFD的验证主要针对计算模型(实现物理概念的计算方法)，诸如网格精度、时间步长以及离散格式等方面的验证。CFD的确认工作主要针对概念模型展开，概念模型即描述物理真实的数理方程和参数，主要包括质量、动量、能量守恒偏微分方程和这些方程的边界条件与初始条件，以及湍流模型等，至于液固流化床内的多相流动，概念模型还包含相间作用力模型(主要是曳力模型)、颗粒相本构关系以及颗粒-壁面作用关系等；确认工作强调的是模型正确求解，一般需要将计算结果和实验数据进行比较。

综上，一个CFD计算结果被认可是需要经过多方面确认的，现实中选择对流动特性产生显著影响的一个或几个概念模型开展确认。具体到液固流化床的CFD模拟，曳力模型的确认工作是必不可少的，这是因为：曳力是最重要的相间作用，实现了动量传递与流体对颗粒的输运；诸多研究表明，使用不恰当的曳力模型会导致液固流化床CFD计算结果与物理真实偏离较大。在确认过程中，往往需要测试若干个曳力模型，挑出计算结果与实验值最为接近的一个，则该结果最为可信。而且，由于床内液固体系的非线性和流域的多态性，当前尚未得到普适的曳力模型，因此，当液固流化床内物性参数和操作条件一旦改变，CFD曳力模型的确认工作就得重新开展。

目前，液固流化床CFD模拟中曳力模型的确认存在以下问题：曳力模型的确认占据整个CFD验证和确认工作相当大的比重，耗费了大量的时间与计算资源，且随着待确认曳力模型数量的增多而极大增多；另一方面，曳力模型始终是多相流学术焦点之一，新模型层出不穷，然而，基于有限资源，通常只对2至3个曳力模型开展确认工作，如此易导致因“样本”数量过少而遗漏了最佳模型。

发明内容

为了解决上述问题，本发明提出了一种液固流化床CFD曳力模型的确认方法。本发明提供的曳力模型确认方法无需开展CFD计算，在确保准确性的前提下，可以节省巨量的计算资源与时间成本，实现即好且快；基于此，该方法可极大地扩充“样本”数量，最终实现对任意数量曳力模型的CFD确认。

为达到上述目的，本发明的具体方案为：

一种液固流化床CFD曳力模型的确认方法，基于对颗粒群体系中单颗粒的受力分析建立力学方程，结合CFD曳力模型表达式，求出床内颗粒体积浓度，并将其与实验测量值进行比较，最终确认出适用于该物性参数和操作条件下的最佳CFD曳力模型。

进一步的，所述方法具备包括如下步骤：

步骤一、给定初始条件；

步骤二、对单个颗粒受力分析，建立力学方程；

步骤三、将步骤一给定的初始条件代入步骤二所得的方程中，解出单颗粒所受曳力；

步骤四、将步骤三所得曳力代入待确认的CFD曳力模型中，解出颗粒体积浓度；

步骤五、将步骤四所得颗粒体积浓度与颗粒体积浓度实验值进行比较，最终确认出适用于该体系下的最佳曳力模型。

进一步的，所述步骤一的初始条件包括：物性参数：颗粒直径d_p，颗粒密度ρ_p，流体密度ρ_f，流体动力粘度μ；操作条件：表观速度u₀；以及在表观速度u₀条件下，实验测得的颗粒体积浓度φ。

进一步的，在所述步骤二中，液固流化床内的单个颗粒仅受重力G、浮力F_b和曳力F_d，其平衡方程为：

G＝F_b+F_d (1)。

进一步的，在所述步骤三中，对于颗粒所受浮力的表述为：颗粒在悬浮液中所受浮力等于其自身体积与悬浮液混合密度的乘积，即：

因此，单个颗粒所受曳力的表达式为：

其中，F_b和F_d分别为浮力和曳力，d_p代表颗粒直径，ρ_p代表颗粒密度，ρ_f代表流体密度，φ为实验所得颗粒体积浓度，g代表重力加速度。

进一步的，在所述步骤四中，曳力模型表征的是单位体积内颗粒受到的总曳力；因此需要先计算单位体积内的颗粒个数N，再将N与步骤3所得的单颗粒曳力相乘；详细计算过程如下：

NF_d＝βu_r (4)

其中，u_r和u₀分别为相间滑移速度和表观速度；β为动量交换系数，是关于颗粒体积浓度和颗粒雷诺数的函数，具体形式由曳力模型给出；将得到的F_d、N和u_r代入至曳力模型，即动量交换系数β，最终解出颗粒体积浓度。

进一步的，包括以下步骤：

1)给定初始条件，其包括：物性参数：颗粒直径d_p，颗粒密度ρ_p，流体密度ρ_f，流体动力粘度μ；操作条件：表观速度u₀；以及在表观速度u₀条件下，实验测得的颗粒体积浓度φ；

2)对单个颗粒进行受力分析，一般将液固流化床床层近似为液固悬浮液，因此床内颗粒仅受重力、浮力和流体所施加的曳力，且处于三力平衡状态；进一步建立力学方程：

G＝F_b+F_d (1)

其中G、F_b和F_d分别表示重力、浮力和曳力；

3)将步骤1所给定的初始条件代入步骤2所建立的力学方程中，浮力的表述为：颗粒在悬浮液中所受浮力等于其自身体积与悬浮液混合密度的乘积：

由此解出单个颗粒所受的曳力为：

其中，F_b和F_d分别为浮力和曳力，d_p代表颗粒直径，ρ_p代表颗粒密度，ρ_f代表流体密度，φ为实验测得的颗粒体积浓度，g代表重力加速度；

4)将步骤3所得曳力代入待确认的CFD曳力模型表达式中，需要特别强调的是，CFD曳力模型表征的是单位体积内颗粒受到的总曳力，因此需要先计算单位体积内的颗粒个数，再将颗粒个数与步骤3所得的单颗粒曳力相乘；其计算过程如下：

NF_d＝βu_r (4)

其中，F_d表示曳力，N代表单位体积内的颗粒个数，u_r和u₀分别表示相间滑移速度和表观速度，d_p代表颗粒直径，φ代表实验测得的颗粒体积浓度；β为动量交换系数，是关于颗粒体积浓度和颗粒雷诺数的函数，具体形式由CFD曳力模型给出；将得到的F_d、N和u_r代入至CFD曳力模型，即动量交换系数β，最终解出颗粒体积浓度；

5)将步骤4中CFD曳力模型计算出的颗粒体积浓度与实验测得的颗粒体积浓度φ进行比较，偏差最小者即为适用于该体系的最佳CFD曳力模型。

进一步的，在步骤四中，以颗粒体积浓度为自变量的CFD曳力模型通常为超越方程，因此，需要采用迭代法来求解颗粒体积浓度。

进一步的，所述曳力模型具体为：Gidaspow曳力模型、Gibilaro曳力模型、DiFelice曳力模型、Syamlal and O’Brien曳力模型以及BVK曳力模型和Rong曳力模型。

进一步的，所述曳力模型具体的表达式如下：Gidaspow曳力模型：

Syamlal and O’Brien曳力模型：

v_r＝0.5(A-0.06Re_p)+0.5[(0.06Re_p)²+0.12Re_p(2B-A)+A²]^0.5 (13)

A＝(1-ε_s)^4.14 (14)

其中Re_p的计算同式(10)；

Gibilaro曳力模型：

Di Felice曳力模型：

其中Re_p的计算同式(17)；

BVK曳力模型：

其中Re_p的计算同式(17)；

Rong曳力模型：

其中Re_p的计算同式(17)；

以上所有曳力模型中，Re_p代表颗粒雷诺数，C_D代表曳力系数，ε_s代表颗粒体积浓度，除ε_s以外，其余变量均由步骤一中所给定的条件确定。需要注意的一点，求解颗粒体积浓度时会遇到超越方程，所以需要采用迭代法来反复逼近真值。

通过采用以上的技术方案，本发明具有以下的有益效果：本发明提供的曳力模型确认方法无需开展CFD计算，在确保准确性的前提下，可以节省巨量的计算资源与时间成本，实现即好且快；在实际应用中效率更高，同时能保证效果，基于此，该方法可极大地扩充“样本”数量，最终实现对任意数量曳力模型的CFD确认。

附图说明

图1为本发明一种确认液固流化床CFD曳力模型方法的流程图。

图2为本发明实施例1的计算结果相对误差图。

图3为为本发明实施例2的计算结果相对误差图。

具体实施方式

下面结合附图对本发明提出的一种液固流化床CFD曳力模型的确认方法进行详细说明。

图1为本发明所公开的一种液固流化床CFD曳力模型的确认方法的流程图，如图1所示，本发明所提出的曳力模型确认方法是应用于液固流化床CFD模拟的，CFD模拟结果需要依靠实验数据来检验可信度，因此取用了相关文献“Chemical Engineering Science 60(2005)1889-1900”和“Chemical Engineering Science 62(2007)6334-6348”的实验数据作为参考。

实施例1

步骤1.给定初始条件：

物性参数：颗粒直径d_p＝0.003m，颗粒密度ρ_p＝2500kg/m³，流体密度ρ_f＝1000kg/m³，流体动力粘度μ＝0.001Pa·s；操作条件：表观速度u₀＝0.07m/s；以及在表观速度u₀＝0.07m/s条件下，实验测得的颗粒体积浓度φ＝0.44。

步骤2.对单个颗粒受力分析，建立力学方程：一般将液固流化床床层近似为液固悬浮液，液固悬浮液有以下两大特点：悬浮液中的颗粒处于力的平衡状态；颗粒间的相互作用可被忽略；因此，可认为液固流化床内的单个颗粒仅受重力G、浮力F_b和流体所施加的曳力F_d，其平衡方程为：

G＝F_b+F_d (1)。

步骤3.对于颗粒所受浮力F_b的表述，采用流态化理论中Gibilaro等人的观点：颗粒在悬浮液中所受浮力等于其自身体积与悬浮液混合密度的乘积，即：

因此，得到单个颗粒所受曳力的表达式：

将步骤1给定的初始条件代入公式(3)，其中，d_p＝0.003m，ρ_p＝2500kg/m³，ρ_f＝1000kg/m³，φ＝0.44，g＝9.81m/s²，解出单颗粒所受曳力F_d＝1.165×10^-4。

步骤4.将步骤3所得曳力代入待确认的CFD曳力模型中，需要特别指出，曳力模型表征的是单位体积内颗粒受到的总曳力，因此需要先计算单位体积内的颗粒个数N，再将N与步骤3所得的单颗粒曳力相乘；详细计算过程如下：

NF_d＝βu_r (4)

最终得到β的表达式为：

其中，u_r和u₀分别为相间滑移速度和表观速度；β为动量交换系数；将由步骤3得到的F_d、以及步骤1中给定的d_p、u₀和φ代入公式(7)，可得到β的值。另一方面，β是关于颗粒体积浓度和颗粒雷诺数的函数，具体表达式由曳力模型给出。

本案例详细对Gidaspow曳力模型、Gibilaro曳力模型、Di Felice曳力模型、Syamlal and O’Brien曳力模型以及BVK曳力模型和Rong曳力模型展开确认，其中，Gidaspow曳力模型、Gibilaro曳力模型、Di Felice曳力模型和Syamlal and O’Brien曳力模型属于较早开发且业已得到广泛应用的曳力模型，而BVK曳力模型和Rong曳力模型属于近期发表的研究成果，尚未在液固流化床的CFD模拟中大量使用。它们的表达式如下：

Gidaspow曳力模型：

Syamlal and O’Brien曳力模型：

v_r＝0.5(A-0.06Re_p)+0.5[(0.06Re_p)²+0.12Re_p(2B-A)+A²]^0.5 (13)

A＝(1-ε_s)^4.14 (14)

其中Re_p的计算同式(10)；

Gibilaro曳力模型：

Di Felice曳力模型：

其中Re_p的计算同式(17)；

BVK曳力模型：

其中Re_p的计算同式(17)；

Rong曳力模型：

其中Re_p的计算同式(17)。

以上所有曳力模型中，Re_p代表颗粒雷诺数，C_D代表曳力系数，ε_s代表颗粒体积浓度，即为所求。观察曳力模型可知，除ε_s以外，其余变量均由步骤1中所给定的条件确定。需要注意的一点，以颗粒体积浓度ε_s为自变量时，曳力模型通常为超越方程，因此需要采用迭代法来反复逼近真值，最终解出颗粒体积浓度ε_s。迭代过程略。

步骤5.将经由步骤4中式(4)—(25)计算所得的颗粒体积浓度ε_s与颗粒体积浓度实验值φ进行比较，最终确认出适用于该物性参数和操作条件下的最佳曳力模型。为了验证本发明所提出方法的有效性，特将本发明所提出的CFD曳力模型确认方法的计算结果与CFD计算结果对比印证，CFD计算由CFD软件实现，须使其所采用的曳力模型和初值条件与本案例相同。计算结果详见表1。

由表1可知：Gibilaro曳力模型得到的颗粒体积浓度最高，Gidaspow曳力模型得到的颗粒体积浓度最低，其余四种曳力模型的结果居于中间。虽然本发明所提出方法得到的颗粒体积浓度与CFD计算结果的绝对数值不同，例如，Gibilaro曳力模型的结果对应本方法和CFD计算分别为0.517和0.473，但其皆大于实验值，Gidaspow曳力模型的结果对应本方法和CFD计算分别为0.368和0.395，但其皆小于实验值，同样地，其余曳力模型的计算结果也呈现这一趋势。

表1案例1不同曳力模型的计算结果对比

为了进一步定量比较各个曳力模型计算结果的准确度，本专利使用相对误差来分析。

其中，ε_s,th代表本专利所计算出的颗粒体积浓度，φ代表颗粒体积浓度的实验值，ε_s,CFD代表CFD计算出的颗粒体积浓度，e_th表示本专利计算结果的相对误差，e_CFD表示CFD计算结果的相对误差。通过计算式(26)和(27)，得到各个曳力模型相对实验值的误差，详见图2。

在图2中，Gb、Gd、S&O和DF分别代表Gibilaro曳力模型、Gidaspow曳力模型、Syamlal&O’Brien曳力模型以及Di Felice曳力模型。观察e_th，BVK曳力模型计算结果的相对误差最小，仅为0.2％，Rong曳力模型稍稍次之，因此，依据本专利所公开的CFD曳力模型确认方法的计算结果来看，适用于该体系下的、推荐使用的曳力模型为BVK曳力模型和Rong曳力模型。观察e_CFD，可以得到相同的结论，如此就直接证明本专利公开的一种液固流化床CFD曳力模型的确认方法是有效的。

实施例2

步骤1.给定初始条件：

物性参数：颗粒直径d_p＝1.13mm，颗粒密度ρ_p＝2540kg/m³，流体密度ρ_f＝999.5kg/m³，流体动力粘度μ＝0.00124Pa·s；操作条件：表观速度u₀＝0.0381m/s；以及在表观速度u₀＝0.0381m/s条件下，实验测得的颗粒体积浓度φ＝0.342。

步骤2-4与案例1实施过程相同，不再赘述。

以下只展示步骤的相关计算结果：由表与图中可知，观察e_th，BVK曳力模型计算结果的相对误差最小，Gidaspow曳力模型稍次之；而观察e_CFD，则有Gidaspow曳力模型相对误差最小，BVK曳力模型稍次之；因此，两种方法得到了相同的结论，即：适用于该体系下的、推荐使用的曳力模型为Gidaspow曳力模型和BVK曳力模型。所以，直接证明本专利公开的一种液固流化床CFD曳力模型的确认方法是有效的。另一方面，对以上6个CFD曳力模型开展CFD的确认计算，在小型服务器上大约耗时4至5天；而使用以上技术方案，在3至4小时内即可完成。因此，通过本发明所公开之方法，可以节省巨量的计算资源与时间成本，实现即好且快。

表2案例2不同曳力模型的计算结果对比

Claims (6)

1.一种液固流化床CFD曳力模型的确认方法，其特征在于，基于对颗粒群体系中单颗粒的受力分析建立力学方程，结合CFD曳力模型表达式，求出床内颗粒体积浓度，并将其与实验测量值进行比较，最终确认出适用于物性参数和操作条件下的最佳CFD曳力模型；

具体包括如下步骤：

步骤一、给定初始条件；

步骤二、对单个颗粒受力分析，建立力学方程；

步骤三、将步骤一给定的初始条件代入步骤二所得的方程中，解出单颗粒所受曳力；

步骤四、将步骤三所得曳力代入待确认的CFD曳力模型中，解出颗粒体积浓度；

步骤五、将步骤四所得颗粒体积浓度与颗粒体积浓度实验值进行比较，最终确认出适用于该体系下的最佳曳力模型；

所述步骤一的初始条件包括：物性参数：颗粒直径d_p，颗粒密度ρ_p，流体密度ρ_f，流体动力粘度μ；操作条件：表观速度u₀；以及在表观速度u₀条件下，实验测得的颗粒体积浓度φ；

所述步骤二中，液固流化床内的单个颗粒仅受重力G、浮力F_b和曳力F_d，其平衡方程为：

G＝F_b+F_d (1)

其中G、F_b和F_d分别表示重力、浮力和曳力。

2.根据权利要求1所述的方法，其特征在于，在所述步骤三中，对于颗粒所受浮力的表述为：颗粒在悬浮液中所受浮力等于其自身体积与悬浮液混合密度的乘积，即：

因此，单个颗粒所受曳力的表达式为：

其中，F_b和F_d分别为浮力和曳力，d_p代表颗粒直径，ρ_p代表颗粒密度，ρf代表流体密度，φ为实验测得的颗粒体积浓度，g代表重力加速度。

3.根据权利要求1所述的方法，其特征在于，在所述步骤四中，曳力模型表征的是单位体积内颗粒受到的总曳力；因此需要先计算单位体积内的颗粒个数N，再将N与步骤3所得的单颗粒曳力相乘；详细计算过程如下：

NF_d＝βu_r (4)

其中，u_r和u₀分别为相间滑移速度和表观速度；β为动量交换系数，是关于颗粒体积浓度和颗粒雷诺数的函数，具体形式由曳力模型给出；将得到的F_d、N和u_r代入至曳力模型，即动量交换系数β，最终解出颗粒体积浓度。

4.根据权利要求3所述的一种液固流化床CFD曳力模型的确认方法，其特征在于：在步骤四中，以颗粒体积浓度为自变量的CFD曳力模型通常为超越方程，因此，需要采用迭代法来求解颗粒体积浓度。

5.根据权利要求3所述的方法，其特征在于：所述曳力模型具体为：Gidaspow曳力模型、Gibilaro曳力模型、Di Felice曳力模型、Syamlal and O’Brien曳力模型以及BVK曳力模型和Rong曳力模型。

6.根据权利要求5所述的方法，其特征在于：所述曳力模型具体的表达式如下：Gidaspow曳力模型：

Syamlal and O’Brien曳力模型：

v_r＝0.5(A-0.06Re_p)+0.5[(0.06Re_p)²+0.12Re_p(2B-A)+A²]^0.5 (13)

A＝(1-ε_s)^4.14 (14)

其中Re_p的计算同式(10)；

Gibilaro曳力模型：

DiFelice曳力模型：

其中Re_p的计算同式(17)；

BVK曳力模型：

其中Re_p的计算同式(17)；

Rong曳力模型：

其中Re_p的计算同式(17)；

以上所有曳力模型中，Re_p代表颗粒雷诺数，C_D代表曳力系数，ε_s代表颗粒体积浓度，除ε_s以外，其余变量均由步骤一中所给定的条件确定。

Images