仿真笔记——Fluent中的多相流模型详解
- 拉格朗日方法——流体被处理为连续相，直接求解时均Navier-Stokes方程；计算流场中大量的粒子，气泡或液滴的运动轨迹，得到离散相的分布规律。离散相和流体相之间可以有动量、质量和能量的交换。基本假设：作为离散的第二相的体积比率很低。
- 欧拉方法——不同的相被处理成互相贯穿的连续介质。引入相体积率的概念，各相的体积率之和等于1。不同的相均满足守恒方程。从实验数据建立一些关系式，使方程组封闭。在Fluent中，有三种欧拉-欧拉多相流模型：流体体积模型（VOF）【分层的或自由表面流】，混合物模型【流动中有相的混合或分离】，欧拉模型【散相的体积分数超过10%】。
VOF模型：通过求解单独的动量方程和处理穿过区域的每一流体的容积比来模拟两种或三种不能混合的流体。典型的应用——流体喷射、流体中大气泡的运动、流体在大坝坝口的流动、气液界面的稳态和瞬态处理等。Mixture模型：一种简化的多相流模型，用于模拟各相有不同速度的多相流，但是假定了在短空间尺度上局部的平衡，相之间的耦合很强。也用于模拟有强烈耦合的各向同性多相流和各相以相同速度运动的多相流。典型的应用——沉降（sedimentation）、气旋分离器、低载荷作用下的多粒子流动、气相容积率很低的泡状流。Eulerian模型：可以模拟多相分离流及相互作用的相（液体、气体、固体），与离散相模型Eulerian-Lagrangian方案只用于离散相不同，在多相流模型中Eulerian方案用于模型中的每一相。材料密度：颗粒在密实状态下，单位体积所具有的质量， 
e=1 —— 弹性碰撞，颗粒碰撞后完全恢复变形，机械能没有损失；e=0 —— 塑性碰撞（完全非弹性碰撞），两颗粒碰撞后不再分开，碰撞引起的变形完全保留下来；0<e<1——实际颗粒的碰撞（非完全弹性碰撞），碰撞过程中有机械能的损耗，对于颗粒运移规律有影响。当量粒径：颗粒形状一般不规则，通常定义一个当量粒径作为颗粒大小的度量，其方法依颗粒大小不同而异。
（圆球体积公式）形状：颗粒整体的几何形态，以球形为标准，定义球度系数来度量颗粒的不同形状。颗粒在运动过程中受到许多种力的作用，不同的力在颗粒运动中起到的作用不同，地位不同，因而处理的方法也不同。4. 压力梯度力，由流场中压力梯度引起的作用力，与惯性力相比，数量级很小，可忽略不计。5. 虚假质量力——附加质量（Added Mass），特例：圆球的附加质量力是惯性力的一半。6. Basset力，发生在粘性流体中，与运动的不稳定性有关8. Saffman升力，流场中存在速度梯度，颗粒受到的升力作用。在速度边界层中，该力的影响比较明显。9. 热泳力，光电泳力，声泳力：在有温度梯度的流场中，使颗粒从高温区向低温区运动的力通常称为热泳力。颗粒吸收光能并加热附近的气体分子，产生类似于热泳力的光电泳力。在声场中，颗粒随着气体振动作用而产生漂移运动。通常情况下，光电泳力和声泳力可忽略不计。10. 颗粒所受的静电力，带有电荷的颗粒在运动中将受到静电力的作用。如果流场中有多个颗粒同时存在，颗粒之间就会发生相互作用。一个颗粒的尾流范围往往比它本身体积大2~3个量级。因此，即使颗粒浓度很低，也存在显著的相互作用——通过流体的间接作用，对颗粒的阻力造成显著影响。作为一个极端的例子，当颗粒一个跟着一个运动时，每个颗粒所受到的阻力比单个颗粒运动受到的阻力小很多。通过实验可以测定各种条件（不同的粒径，不同颗粒浓度等）下的表观阻力系数。例如，用照相或其它方法测出某些颗粒的平均速度及其变化后，就可以求出它们的雷诺数Re和阻力系数CD。很多作者用实验方法研究过颗粒群的阻力，并归纳出计算阻力的经验公式。Rudinger在激波管中用纹影照相和光散射法测量颗粒运动得到 实验本身的误差，颗粒群阻力的测量远比单颗粒阻力测量困难，但这些实验的总的趋势应该说是可信的。 实验条件不同，如流动的湍流度、颗粒粗糙度、以及粒度非均匀性等，颗粒的静电效应、旋转效应、流体的入口条件等，也是造成阻力公式不尽相同的原因。- 气-液相界面的尺寸很小，而且具有与流体一样的、良好的流动性；
- 气-液相界面上存在表面张力，在部分情况下，这种表面张力对两相流的流动状态具有决定性作用；
- 由于相界面的流动和分布特性的不同，气-液两相流具有多种宏观特性完全不同的流型，而泡状流仅仅是气-液两相流的众多复杂流型中的比较简单的一种。
气-液两相流的数值模拟比气-固两相流复杂得多，其流动结构和宏观特性都与气-液相界面的分布有关。按照描述相界面的方法不同，气-液两相流的数值模拟方法包括下述几类。在研究江河等自由表面问题时，最容易想到的描述相界面位置的简单方法就是标高法，即高度函数法——将相界面到一个参考平面（如河底水平面）或参考直线的距离定义为参考面（线）上各点位置的函数，由此得出相界面的分布形状。优点：相界面形状及分布状况简单、明了，尤其适合于描述像明渠流动这样的两相流问题。缺点：当界面斜率超过网格的高宽比 
时，这种方法不能很好的工作；尤其当界面发生波浪式翻卷或破碎时，高度函数变为多值函数，该方法无能为力。PIC（Particle-In-Cell）方法最早由Los Alamos实验室的Harlow和Welch等人提出。 连续流场被看成是由有限个分布在Euler网格内、质量集中在网格中心的流体质点构成的体系，每个流体质点具有各自的质量、动量和能量。各个流体质点完全按照自己的“个性”去运动，（Lagrangian步）对质点的过分“聚集”或“稀疏”的情况进行判断和重新调整，使流体质点在空间中的位置分布比较合理，每一个空间位置必须有，而且只能有一个流体质点对流过网格边界的质量及其携带的动量和能量进行修正，计算出各网格在下一步的速度和内能等。（Euler步）优点：为气-液两相流的数值模拟提供了一种有效途径，为其它方法的产生奠定了基础。缺点：需要对每一个流体质点进行跟踪，并存储每一个空间位置上流体质点的各种运动学和动力学参数，计算工作量、存储空间要求很高。MAC（Marker-and-Cell）方法实际是早期PIC方法的改进，由Harlow和Welch于1965年提出（Harlow & Welch，1965；Daly，1969）。在MAC方法中，如果一个计算单元中含有记号粒子，而其相邻单元不含记号粒子，则这个计算单元在相界面上（含有相界面）；相界面在该计算单元内的具体 位置需根据该单元内记号粒子的分布状况计算确定。确定相界面的这种思想也被后来的VOF方法吸收并得到发展。最大贡献：排除了与界面相交或重叠有关的逻辑判断问题线段法是表示相界面的一种直观方法，即用一连串的线段或一批由线段连接起来的点表示一个相界面（刘全，2002；Lebaigue，1998）。需要存储上述线段的每一个点的坐标，为了保持精确性，相邻两点之间的距离一般不超过最小网格尺度（Dx或Dy），Lagrange方法。缺点：计算过程中需要根据界面的变化情况在部分区域增加线段（点），或在部分区域删除部分线段（点），因而需要做大量的逻辑判断工作，实施过程比较复杂，尤其是相界面发生交叉或折叠时。边界积分法的基本思路：在特定条件下，比如界面变形比较小时，可认为界面附近的液体运动是无旋的，液体运动简化为理想不可压缩流体的无旋运动，其速度势函数F满足Laplace方程。在一定的初始条件和边界条件下，对气-液界面做Lagrangian积分，可得出自由表面随时间变化的控制方程，进一步求解自由面的控制方程，可得出自由界面随时间的变化规律。
边界积分法将气-液自由面用伪弧长参数和样条曲线表示，将流体内部点的速度与界面上的速度和应力联系起来，特别适合于研究二维结构中变化不太大的相界面问题——自由表面的流动与变形问题（罗朝霞，2002；Davidson，2000）。
船舶在波浪中的运动，海洋工程问题。
（1）构造伪弧长参数样条曲线：将自由表面采用伪弧长参数样条曲线来表示，把自由表面分成几段，节点之间的弧长用直线距离来代替，然后用三次样条差值求出相关的函数。（2）用迭代法求解边界积分方程：在迭代求解过程中，一般采用亚松弛，例如松弛因子可取0.8。→用势流理论得到离散方程组。（3）重分网格：在追踪自由面位置的过程中，为了避免网格畸形，必须采用重分网格技术。主要优点：利用势函数将所研究问题的空间维数降低一阶，但又不影响解的精度；能直接计算出界面速度，追踪界面变形，并将表面张力及其它一些表面效应直接包含在计算过程中。局限性：对变系数、非线性问题的适用性不好；计算过程中需要随时对界面进行网格加密或拆分；数值计算方面，积分核的奇异性和由控制方程的离散化得到的代数方程组的非稀疏性，给计算增加了困难；当相界面发生交叉或折叠时，该方法变得复杂。该方法最初由Osher和Sethian（1988，1996）提出，主要应用于智能控制、图像处理、材料微细加工等领域。经过Sethian，Sussman等人的进一步发展，该方法被用于气-液两相流的数值模拟研究。基本思路：把气-液相界面的传播用一个高阶函数f的零点值来表示，由（Level Set函数）f的代数值来区分计算区域中的各相。这里，Level Set函数取代了VOF方法中的相函数F。Level Set函数的特点：该函数始终为距离函数，可以方便的计算相界面的曲率、法向向量等几何参数，从而将相界面上的表面张力用连续函数的形式表示出来。优点：求解思路比较容易理解，相界面可以被表示为连续函数，便于做数学运算；可求解相界面的几何特性参数，从而求解表面张力；该方法也容易向高阶空间推广。缺点：某些情况下，容易造成算法的不收敛或数值不稳定性。1981年，Hirt和Nichols首先用VOF方法对溃坝及Rayleigh-Taylor不稳定性现象进行了模拟，对运动相界面问题的数值模拟研究做出了开创性贡献（Hirt & Nichols，1981；Guegffier，1999）。VOF方法用相函数（Phase Function）F取代了MAC方法中的虚拟无质量粒子，可看作是MAC方法的一个改进，类似于气-液两相流中的截面含气率（或容积含气率），表示某一相介质占据网格面积（二维）或体积（三维）的分数。
计算量大，计算结果的稳定性。目前较为流行的方法：Level Set，VOF。大涡模拟（Large Eddy Simulation，LES）是为了深入研究单相流体的湍流运动而发展出的一种新的数值模拟方法。（1）大尺度旋涡主要影响平均流动，各种变量的湍流扩散，热量、质量及动量的交换以及雷诺应力的产生都是通过大尺度旋涡实现的。（2）小尺度旋涡主要对耗散起作用，通过耗散脉动来影响各种变量。（2）小旋涡可通过亚网格尺度模型建立与大尺度旋涡的关系。目前主要用于模拟由连续流体相与弥散颗粒相组成的两相流动系统。
直接数值模拟（Direct Numerical Simulation，DNS）是在湍流耗散尺度相当的网格上直接求解瞬态的三维Navier-Stokes方程，不需要封闭模型，但要求高精度的数值方法和合适的边界条件。对于两相流动，对有限尺寸和有限数目的颗粒/气泡的直接模拟，可以直接得到颗粒/气泡在流体中的受力及其尾涡对流体湍流的作用。直接模拟的计算量很大，计算时间很长，一般需要在大型计算机上完成，目前主要用于低Reynolds数和简单小尺寸流动的模拟，尚难以直接用于工程问题。（10）多相流的Lattice-Boltzmann方法前面提到的多相流的数值模拟方法都是在连续介质假设的前提下展开的，需要先根据连续介质假定，建立起流体运动所遵循的微分方程，然后，以此微分方程为出发点，采用有限差分、有限体积、有限元或有限谱等离散方法，对微分方程进行离散，再用适当的数值方法进行求解。分子动力学类方法根据分子运动理论直接建立起简化的动力学模型，将流体的宏观运动看作是大量流体分子微观运动的统计平均的结果。介于连续介质模型和分子动力学模型之间的格子-粒子类方法：将流体的运动场用一系列规则的格子表示，计算就在这些格子（Lattice）上进行。格子-粒子类方法：将流体的运动场用一系列规则的格子表示，计算就在这些格子（Lattice）上进行。这些格子的尺度远比分子的平均自由程大，但又比有限差分的步长或有限体积法中的控制体积宽度小许多；
计算对象不再是数目庞大的流体分子个体，而是数目已大大减少的流体“粒子”。这些粒子比分子级别要大，但其质量又比有限体积法中的控制体积内的流体质量要小得多。大大降低了对计算机存储量的要求。
“粒子”只能在格子之间按一定的规则、沿格线以相同的离散时间步长运动，所以“粒子”的运动轨迹也大为简化。
格子-Boltzmann方法：是Boltzmann方程的一种特殊的离散形式，是上述格子-粒子类方法的改进和发展。该方法借用分布函数 
来描述由N个粒子构成的离散系统。这里的qi和pi分别代表第i个粒子的广义坐标和动量，分布函数fN则包含了所有动力学过程的全部统计信息，Boltzmann方程则描述流体系统的粒子分布函数的演化过程。通过积分，可以从分布函数得到流体的宏观密度、动量和能量。该方法完全不同于传统数值方法，它的微观粒子背景使得它具有许多其他数值方法所没有的独特优点。除了上述方法，离散涡模拟也是一种重要的多相流数值计算方法，离散涡（Discrete Vortex）又叫随机涡（Random Vortex）模拟，是把湍流流场分成一系列大尺度的涡元，用涡元的随机运动来模拟湍流。横向冲刷钝体的流体在一定条件下会产生旋涡脱落现象，在尾流中形成所谓的卡门涡街。涡街会使物体表面的压力幅值发生周期性变化，导致物体在与来流垂直的方向上受到大小和方向周期性变化的流体作用力，从而产生所谓的流致振动现象。高空假设的管道、输电系统中过江等架设的大跨度电缆、潜艇的潜望镜、高耸的烟囱、大型冷却塔、天线、电视塔、水塔和桥梁等。如果流致振动的频率与结构的固有频率接近，就会激起结构物的共振，导致严重的破坏。现有研究成果认为，流动引起振动的方式有四种：旋涡分离【比较普遍】、湍流抖振、声共鸣和流体弹性不稳定性。圆柱绕流：对于不可压缩流体，光滑圆柱体的旋涡脱落形式随Re数变化。流动无分离；
尾流中对称的稳定旋涡；
层流涡街；
湍流涡街（亚临界区）；
超临界区；
高超临界区（周期性旋涡脱落和湍流涡街重新建立）。
（2）管束绕流
管束常见的排列方式：
单行管，正方形排列，旋转正方形排列，正三角形排列，旋转正三角形排列。相邻两管中心线之间的距离P与管子外直径d之比，称为节距比（P/d）。
对于两圆柱前后排列的结构，根据S1/d和S2/d的不同，会出现如下情况：(a)两圆柱后方出现单涡街；(b)前一圆柱脱落边界层交替再贴附在后圆柱上；(c)准稳定的边界层再贴附；(d)间歇性旋涡脱落；(e)间断性的从形式d跳到f或反之；(f)双涡街等形式。
对气液两相流体横掠单根柱体旋涡脱落特性的研究无论在实验还是在理论上都研究得很少。- 气液两相流绕流柱体旋涡脱落引起的柱体表面周向压力分布特性等。
（1）单相流动中，Strouhal数是一个表明旋涡脱落特性的相似准则数，用符号St表示，
（2）单相流动中，管束的Strouhal数仍用上式计算，其中，迎流面宽度采用管子直径d，特征速度u的选择：
其中，S1为横向节距。
（3）气液两相流体，横掠柱体时，Strouhal数可采用下式表达
式中，f——气液两相涡街频率（Hz）；
W——涡街发生体特征尺寸（本书中该尺寸为发生体的迎流面宽度）（m）；
um——气液两相混合物流速。
实验条件：两相流型为细泡状，实验段管子横截面，柱体剖面形状、尺寸，流动方向（垂直上升、垂直下降），截面含气率等。当气液两相流横掠单根柱体时，随着旋涡在柱体后部两侧交替脱落，流体作用在柱体表面的静压也随之发生脉动。在细泡状流型的条件下，观察气液两相流横掠管束时的旋涡脱落。（串列双圆柱实验段、并列双圆柱实验段）
气液两相混合物沿水平方向横向冲刷水平布置的错列管束，可划分为四种流型：（a）泡状流（bubbly flow）；（b）分层流（Stratified flow）；（c）分层-喷雾流（Stratified-spray flow）；（d）喷雾流（spray flow）。直圆管中垂直下降气液两相流的流型如下。流型研究是对一个气液两相流系统开展研究的起点。
目前，文献中已有大量的有关压力降的模型和关联式，这说明压力降的重要性，但对于两相压降，现在还没有一个特别准确的关联式。原因一：任何一个关联式都不可能考虑到所有的影响因素，例如，进口条件会对压力降产生很大影响，会持续影响到下游几百个直径。原因二：任何一个关联式不可能考虑到所有的流型，流型的不同意味着不同的流体之间的相互作用，从而影响着压力降的大小。两相流压力降计算的基础：质量、动量和能量守恒方程。传统的处理方法，以现象学为基础，得到一组包括下列假设的简化方程：（1）管道的任何流通截面上，流体的压力不变；（2）尽管相间存在速度差，但是流通管道中每一相的速度不变。
当两相流体流过突扩、突缩、孔板、弯头和阀门等部件时，静压变化、能量损失目前还研究得比较少。两相流体在横截面积发生突然增大的管内流动时，将产生流体分离现象，在单相流中所使用的方法，仍然适用于一维的两相分离流动。用于单相流体流过突缩接头时的计算方法，也适用于两相分离流流过突缩接头时的计算。
汽液两相流过锐边孔板时会产生较大的压力降，可以采用常规的单相不可压缩流动的方法来处理分离两相流动。
1. 多相流及其应用，车得福 李会雄 编著，西安交通大学出版社，2007年11月。2. 液-固两相流基础，岳湘安 著，石油工业出版社，1996年4月。3. Fluent培训材料，Fluent 6.1 User’s Guide，Fluent Inc.，20034. 计算流体动力学分析——CFD软件原理与应用，王福军 编著，清华大学出版社，2004年9月。
