（检测技术与自动化装置专业论文）油水两相流流型检测方法研究.pdf

摘要 随着采油技术的发展，油井产液剖面测井也变得日益重要，而流量、含水率等参数 的准确测量往往依赖于对流型的了解；在油气储运过程中，由于管路结垢、腐蚀而造成 的泄漏事故时有发生，这不仅造成了重大的环境污染，而且造成了巨大的经济损失，而 管道的结垢、腐蚀的成因也与流型有很大的关系。因此，油水两相流流型特征的研究具 有重要的实际应用价值。 本文开发了基于格子b o l t z m a n n 理论的油水两相流流型模拟软件，设计了探针测量 油水两相流测量系统，并将两者结合检测了油水两相流的流型，取得了如下研究成果： 1 、本文设计了垂直上升油水两相流测量系统，由探针测量系统获得了三种流型的 探针信号，并分析了垂直油井中三种流型的特点，通过实验发现：过渡流型中油水的分 布具有随机性，油包水流型和水包油流型不仅随着流速的不同而变化，而且随着管道的 深度不同而变化。 2 、通过软件模拟获得了水包油、过渡流型和油包水三种流型的典型模拟图，根据 模拟结果获得了油水两相流的流型转换图。经分析发现：由软件模拟的流型特点和由探 针获得的流型特点具有较好的一致性。 3 、通过对探针法测量信号分析获得的流型转换图和软件模拟的流型转换图进行比 较可知：在水相表观速度较高的情况下，两个流型转换图基本一致，其他区域稍有偏差。 这为以后水平油水油井流型的研究提供了一条路径，并为油、水、气三相流的模拟打下 了坚实的基础。 关键词：油水两相流探针格子b l o t z m a n n 流型识别 a b s t r a c t w i t l lt h ed e v e l o p m e n to ft h eo i lw e l lp r o d u c t i o nt e c h n o l o g y , t h ep r o d u c i n gp r o f i l ew e l l l o g g i n gi sb e c o m i n gm o r ea n dm o r ei m p o r t a n t t h ee x a c t l ym e a s u r i n go ft h ep a r a m e t e ro f f l o w r a t ea n dw a t e rv o l u m ef a r c t i o na l w a y sd e p e n d so nt h eu n d e r s t a n d i n go ft h ef l o wp a t t e r n t h el e a ko f t e nh a p p e n si nt h eo i l & g a ss t o r a g ea n dt r a n s p o r t a t i o nb e c a u s eo fc o r r o s i o na n d s c a l ef o r m a t i o no f p i p e s i tl e a d st on o to n l ye n v i r o n m e n t a lc o n t a m i n a t i o nb u ta l s oe c o n o m i c l o s s t h er e a s o no fc o r r o s i o na n ds c a l ef o r m a t i o no fp i p e sh a ss o m e t h i n gt od ow i t hf l o w p a t t e r n s ot h er e s e a r c ho no i l w a t e rt w o - p h a s ef l o wp a t t e r nh a sas i g n i f i c a n ta p p l i c a t i o n v a l u e i no r d e rt os i m u l a t eo i l w a t e rt w o p h a s ef l o wp a t t e m ，t h es o f t w a r eh a sb e e nd e v e l o p e d p r o b e sm e a s u r e m e n ts y s t e mh a sb e e nd e s i g n e di no r d e rt ov e r i f yt h es i m u l a t i o no fo i l w a t e r t w o - p h a s ef l o wp a t t e r n b yc o m b i n i n gt h ep r o b e sm e a s u r e m e n ts y s t e m 谢t 1 1s i m u l a t i o n ， f o l l o w i n gc o n c l u s i o n sh a v eb e e na c h i e v e d ： 1 u s i n gp r o b e sm e a s u r e m e n ts y s t e m ，t h es i g n a l so ft h r e ek i n d so ff l o wp a t t e mh a v e b e e ng o t t h ec h a r a c t e r i s t i co ft h r e ek i n d so ff l o wp a t t e mo fv e r t i c a lo i lw e l lh a sb e e nf o u n d a c c o r d i n ge x p e r i m e n tr e s u l t s ，w eh a v ed i s c o v e r e dt h a tt h ed i s t r i b u t i o no fo i la n dw a t e ri n o i l - i n w a t e rf l o wi sr a n d o mb u ti ti sd i f f e r e n tw i t hd i f f e r e n td e p t ho rd i f f e r e n tv e l o c i t yi n o i l i n w a t e rf l o wa n dw a t e r - i n o i lf l o w 2 t h r e ek i n d so ff l o wp a t t e r n ，w h i c ha r eo i l - i n - w a t e rf l o w , t r a n s i t i o n a lf l o wa n d w a t e r - i n o i lf l o w , h a v eb e e ng o tb ys i m u l a t i o n a c c o r d i n gt ot h er e s u l to fs i m u l a t i o n ，t h e o i l w a t e rt w o - p h a s ef l o wp a t t e mt r a n s i t i o nb o u n d a r ym o d e lh a sb e e ng o tb y b yt h ea n a l y s i s o fs i m u l a t i o n ，t h ec h a r a c t e r i s t i co ft h r e ek i n d so ff l o wp a t t e r no fv e r t i c a lo i lw h i c hh a sb e e n g o tb ya n a l y s i so ft h es i g n a l si sc o n s i s t e n tw i t hr e s u l t sb ys i m u l a t i o n 3 t h et w ov e r t i c a lu p w a r do i l - w a t e rt w o p h a s ef l o wp a t t e mt r a n s i t i o nb o u n d a r ym o d e l s a r ec o n s i s t e n tw h e nw a t e rs u p e r f i c i a lv e l o c i t yi sh i g h ，w h i l et h e r ei sd e v i a t i o nw h e nw a t e r s u p e r f i c i a lv e l o c i t yi sl o w i tp r o v i d e san e ww a yf o rt h er e s e a r c ho fh o r i z o n t a lo i l - w a t e r t w o p h a s ef l o w , a n di ti st h ef o u n d a t i o no fo i l w a t e r - g a st h r e e p h a s er e s e a r c h k e yw o r d s ：o i l - w a t e rt w o - p h a s ef l o w , p r o b e ，l a t t i c eb l o t z m a n n ，f l o wp a t t e r ni d e n t i f i c a t i o n 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取 得的研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他 人已经发表或撰写过的研究成果，也不包含为获得 墨盗墨兰太至 或 其他教育机构的学位或证书而使用过的材料。与我二同工作的同志对本研 究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。 学位论文作者签名：韩幽岛 签字日期：沁7 年 月y 日 学位论文版权使用授权书 本学位论文作者完全了解 墨盗墨兰盘茎有关保留、使用学位论文 的规定。特授权墨盗墨兰太堂 可以将学位论文的全部或部分内容编入 有关数据库进行检索，并采用影印、缩印或扫描等复制手段保存、汇编， 以供查阅和借阅。同意学校向国家有关部门或机构送交论文的复本和电子 文件。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：韩也结 签字日期：枷? 年弓月日 导一：教辩 特啪：7 一日 第一章绪论 第一章绪论 1 1 油水两相流流型检测技术研究意义 相的概念通常是指系统中具有相同成份及相同理化性质的均匀物质部分。当流体中 同时存在不同相成份时就形成了多相流。如同时存在气相和液相两种流动介质的气液两 相流；同时存在气相和固相的气固两相流，以及同时存在油和水两种互不相容液相组分 的油水两相流等。 随着现代工业的发展及两相流在工业中应用越来越广泛，对生产过程中两相流的计 量、控制要求越来越高，两相流参数的准确测量变得至关重要。揭示两相流流动特性、 建立两相流模型并为预测流动过程、设计和优化控制提供必要的依据，对涉及两相流的 大型系统安全高效运行具有指导意义。尤其在石油勘探开发及油气运输等领域，普遍存 在着油水两相流动问题【l 】。我国在“十五”科研项目及“8 6 3 等高科技项目中都对油水 两相流研究给予了高度重视，并在海上及陆地定向井及近水平井产液剖面测试技术中取 得较大进展，但由于油水两相流流动形态的多态性及复杂性，致使其流动参数检测水平 远未满足石油工业应用要求。 在石油开发领域，油井多相流测井资料解释精度和油水两相流动特性密切相关，它对 油井生产特性及动态监测影响至关重要。同时，随着油气田开发的不断深入，我国陆上 大多数注水开发油田已经进入了中高含水期。这时期，油气储集层和油井内的油、水 分布发生很大变化，此时的油水流型比低含水期的油水流型要复杂的多。为了掌握油气 储集层及油井内流动特性变化，油田动态监测系统要求及时地了解各个油气生产层的状 态，准确地确定油井产出剖面，以便合理地调整油气开发方案，使油井处于最佳或正常 的生产状态，最终达到提高油田采收率目的。 由于油井多相流系统的流动特性十分复杂，单从流动机理来分析和解决油井产出剖 面问题是十分困难的。为了实现上述储集层的动态监测需要，必须深入研究油井多相流 系统，建立系统模型，研究先进的油井多相流参数测量技术，对优化油井生产动态特性 及正确地指导油田开发都具有极其重要的意义。 为了实现控制和预测油水两相流动，除了要知道油水两相流流动条件、流体性质和 流体组分外，流型特征的识别，也是其必不可少的一个重要条件。在两相流的流动过程 中，由于流动介质除与管壁之间存在着相互作用外，两相界面之间也存在着互相作用。 这不但影响两相流动的力学关系，而且影响传热和传质性能，两相流动中各相界面间的 相对速度、相界面间的扰动等使流型多种多样、十分复杂1 2 j ，因而流型对两相流特征参 数的准确测量产生很大影响，特别是在过渡区产生很大的误差。另外，由于不同流型下 测量仪器的响应差异很大，为了使测试数据能较正确地反映实际流动状况，需要了解井 下油水两相流流型信息，所以，油水两相流流动结构及其相态逆转特性研究是亟待要解 第一章绪论 决问题1 3 1 。 1 2 油水两相流的研究现状 国内外许多学者都对垂直流动进行了研究。油水两相流由于两相介质均为液体，且 密度相近，所以其流动特性与气液两相流有很大的不同，不能按照气液两相流流型分类 方法进行划分i 训。g o v i e r 等【5 j 在内径为2 6 4n u n ，高度为1 1 3m 的垂直透明管中进行了油 水两相流实验，认为液液体系中观察到的流型次序和在气液体系中观察到的相同，两种 体系中流型明显区别是由于粘度、密度和界面张力不同所造成的，而且认为在任意的水 流量下，只要空气的流量大约1 0 倍于油流量时，气一水体系即表现出与油水体系相同的 流型。所给出的油水流型图与气液类似，分为：( 1 ) 泡状流( 油为分散相，水为连续相) ； ( 2 ) 段塞流( 油为分散相，水为连续相) ；( 3 ) 泡沫流( 油泡或水滴分散，无固定连续相) ； ( 4 ) 雾状流( 水为分散相，油为连续相) 。h a s s o n 等人1 6 j 做了等密度油水混合物的实验， 发现壁面的湿润效果对环状流的破裂及从环状流向其他流型转变起着重要的作用。 v i g n e a u x 等 7 1 在内径为2 0c m 垂直上升管中观察到在油水两相流过渡流型转换时其持水 率变化范围为0 2 0 0 3 0 。在油水垂直油水两相流的流型模型中，较为经典的为f l o r e s 的 模型“j ，f l o r e s 利用内径为5 0 8n u t i ，长为3 8 1r n n l 的管路，采用矿物油和水( 油水密度 比o 8 5 ，粘度l 卜, 2 0 0 ，5 0 ) 对研究了垂直上升管中油水两相流流动特性，在垂直管中 油水两相流动可划分为油为连续相和水为连续相两大类( 图1 - 1 ) 。水为连续相流型，包 括油滴分散流型( do w ) 、细小油滴均匀分散流型( v f do w ) 及油相在水中混状流型 ( o wc f ) ；油为连续相，包括水滴分散流型( dw o ) 、细小水滴均匀分散流型( v f dw o ) 和水相在油中混状流型( w oc f ) ，同时提出了流型转变机理模型，指出了在评价细小分 散泡流转变时是基于湍流动能及液滴表面自由能两个概念，而在评价泡沫流及包含相转 变时是基于聚结或聚集的概念。钟兴福等【8 j 在模拟井实验基础上，通过直接观察测量方 法，在1 2 5n l i i l 管径中对垂直上升管中油水两相流流型进行了重新划分，将流型划分为 油泡流、分散油泡流、混状流和水泡流，并且给出了油水两相流流型转化边界，建立了 用流量法和持水率法判别流型标准，实验发现小直径管中的流型划分结果不能成比例地 应用到大直径管中，而且划分标准或界限也不能按比例放大或缩小。赵东建等【9 】应用双 头电导探针测量系统，对垂直上升管内油水两相分散流局部相分布特性进行了系统测 量，得到了油水两相分散流的局部含油率分布类型图，研究结果表明低折算水速和低折 算油速条件下，局部含油率在实验段截面上呈抛物线型局部分布特征，局部最大值出现 在实验段中心区域，随折算油速增大，油滴受到横向力作用，逐渐向实验段壁面区域迁 移，形成局部含油率的壁面峰值分布特性，当折算水速大于0 8m s 时，局部含油率在实 验段截面上呈均匀分布。金宁德等【l o j 基于运动波理论建立了垂直上升管中油水两相流运 动波传播方程，计算得到了管内无仪器插入体情况下油水两相流运动波传播速度特性曲 线，根据曲线特征确定了油水两相流过渡流型存在不同的持水率，将流型分为三种：( 1 ) 水包油流型，油泡分散在水中。水为连续相，持水率大于6 0 ；( 2 ) 过渡流，油水混合在 一起，没有连续相，持水率在2 5 一6 0 之间；( 3 ) 油包水流型，水滴分散在油中，油携 第一章绪论 带水滴流动。油为连续相，持水率小于2 5 。 圄图圉囤 圈l - 1垂直上升管中油水两相流流型 1 3 目前油水两相流的研究存在的问题 虽然对流型的研究多年来一直在进行，但是由于流动条件变化的多样性以及实验与 理论研究角度不同，基于流动结构形态学上的流型划分方法始终没有特别统一的结论。 目前油水两相流流型研究主要存在以下几个问题： l 、流型定义模糊性：目前对于各种流型的定义只是建立在主观观察的结果上，不 可避免的存在观察者主观因素的干扰。 2 、流型定义多样性：对于各种流型，目前还没有建立起一套统一的公认定义。 3 、流型定义定性特性：对于流型判别只能定性判断，公认的定量判断方法还没有。 1 4 油水两相流测量方法现状 油水两相流的流型1 份复杂，在液液两相流中分别存在液一液两相分界面，由于相 问界面形状都可以发生形变，因此油水两相介质的分布状态随分相流量变化也不断变 化，油水两相流的流型不仅是多种多样的而且其变化常带有随机性。 到目前为止，在工程上确定流型主要是依靠流型图( f l o wp a t t e r nm a p ) 。流型图是通 过流型实验或通过流型计算所取得的流型及其转变与各种参数的坐标关系。流型图上的 区域可以分别表示相应的流型，图中的分界线，是两种流型的转变区。把这些分界线用 流型图中的坐标参数拟和，则得到描述分界线的函数式，也称为流型判别式。 除去依靠流型图可以取得流型判别式外，还可以根据流型形成的物理条件、几何条 件和力学条件，进行组合并关联一些无园次准则，从而得到更复杂的流型判别式。现有 流型图的坐标参数都比较简单，它不可能充分描述与两相流众多参数有关的流型。因此， 由无园次准则取得的判别式，往往更为精确，一般可直接用于工程计算程序的编排中， 而流型图则需预先经过拟合方能纳入程序。 根据流型图和流型转换判据虽然可以确定流型，但最有效的方法还是利用某种仪表 第一章绪论 直接测量管道内两相流的实际流型。根据工作原理，测量方法可以分为直接测量法和间 接测量法两类【l 。其中直接测量法就是根据两相流的流动形式，直接确定流型，如目测 法、高速摄影法、x 射线法以及探针法等。间接测量法是通过对反映两相流波动特性的 信号进行处理分析，来提取流动状况的基本特征，进而确定流型。下面我们主要讨论用 间接测量法进行流型识别的发展情况。 l 、电导测量方法 这是一种以两相流或多相流为主要对象的在线实时检测技术。它能进行两相流相含 率和流体流动速度的测量，并可经过进一步提取处理若干两相流的特征参数。其主要测 量方法包括纵向多极阵列电导法【1 2 , 1 3 j 。电导测量方法在油井测量中已经投入实际使用。 2 、幅值域处理方法 幅值域的处理指对信号幅度进行统计处理，所给出的结果可以反应信号幅度变化的 分布规律，常用的参数有概率密度函数、方差、均值等。由于气液两相流参数直接影响 着流动信号的统计特性，对于测量获得的气液两相流体系中的常规信号，比如局部空隙 率、压力、压差信号及电导、电容信号在幅值域上进行统计处理往往可以实现对流型、 空隙率等流动参数的检测估计。 j o n e s 和z u b e r l l 4 】对x 射线得到的空隙率信号进行了概率密度函数的幅值分析，并 用于辨识油水两相流的流型。v i n c e 和l a h a y1 1 5 j 进一步用p d f 的方差来进行流型判别。 s o n g 等【l6 】发现泡状流到段塞流的转化可以通过p d f 函数的形状得到很好的反映。应用 概率密度函数进行流型辨识的方法主要是依靠对p d f 函数的峰值形状的判别。但实际 应用中，许多时候p d f 函数的形状并不是很典型，不能很容易的确定峰值的特征，这 就给流型辨识工作带来了一定的困难。而且，为了保证p d f 函数数据的准确性和正确 性，数据长度是越长越好，对于非平稳的随机过程更是如此，这就不利于流型辨识的实 时性要求。 3 、频域处理方法 早在1 9 6 6 年h u b b a r d 和d u k l e r t l 。7 】就将水平管道气液两相流压力信号的功率谱分析 结果用于流型判别。功率谱密度函数分析法的一个关键性问题是功率谱密度的分布不完 全取决于流型，而且与两相流的流速关系较大，但流速一般难以预知，因此这种分析法 未z 月匕e , 一z t k e t 好的在工程中得到应用。而且使用功率谱密度函数分析法的一个隐含的假设是， 所分析的信号必须是平稳的，而平稳性假设在两相流对象中往往不能得到满足，这必然 为识别结果带来误差。但是一些研究结果表明【1 8 1 9 】，这种分析方法在预估流型转换方面 还是十分有效的。 4 、软测量方法 软测量技术是依据某种最优化准则，利用过程中易测量的辅助变量( a u x i l i a r y v a r i a b l e ) ，通过构造某种数学关系和估算方法，实现对难以直接测量的主导变量( p r i m a r y v a r i a b l e ) 的在线估计。 在多相流研究领域，目前常用的软测量方法主要有以下几种【2 0 】：( 1 ) 基于机理分析的 软测量方法；( 2 ) 基于人工智能的软测量方法；( 3 ) 基于过程层析成像的软测量方法；( 4 ) 基于动态信号分析的软测量方法。 第一章绪论 软测量方法已经在油水两相流领域得到了很好的应用，吴浩江等【2 l 】采用人工神经网 络对油水两相流流型进行了辨识。赵鑫等 2 2 j 对垂直上升管中油水两相流流动测得的电导 波动信号，从时域和频域内提取了1 1 个反映油水两相流流动特性的特征量。将这些特征 量作为人工神经网络的输入量，在总流量1 0 6 0m 3 d a y 及含水率5 l 一9 1 范围内，采用 软测量方法较好地实现了油水两相流含水率预测。 5 、非线性分析方法 近年来随着数学理论的发展，非线性理论越来越为人们所重视。混沌、分形等理论 正是分析和处理非线性系统的有效数学工具。 油水两相流系统从自身特点分析，是一个非线性耗散动力学系统，它呈现出一类混 沌行为，而且可以用分形特性来描述多相流系统的流型，从而产生了混沌和分形 ( c h a o s & f r a c t a l ) 分析流型的方法，它以全新的角度认识流型的特征，被认为是比较有发 展前景的一条途径。 进入二十世纪九十年代，应用分形、混沌理论来研究多相流的成果日趋增多。f a n 等1 2 3 ，2 4 j 利用差压波动信号研究了气一液一固及液一固流化床，研究结果表明，从时间序 列提取的分形维数与流化床流动特性有较好的对应关系。d a w 等【2 5 】研究了复杂气固两相 流的混沌特性，通过对差压波动信号的混沌分析，得出了复杂气固两相流具有低维的奇 异吸引子。对于垂直管中的油水两相流，o d d i e 等【2 6 l 利用电容、电导和超声技术研究了 混沌吸引子相关维数均近似为5 。金宁德等【2 7 郐】利用符号统计学及复杂性测度理论，通 过非线性时间序列分析来寻求流型辨识的客观指示器。 采用信号处理的方法辨识两相流流型，避免了复杂的机理分析，但特征量的提取方 法仍存在很多问题，一种信号处理方法不能解决不同测量对象及不同工况条件下的所有 流型辨识问题；一些非线性特征量计算量大，提取并非容易，这也限制了该方法的实时 性和推广应用；此外一些算法本身的鲁棒性及全局收敛性仍有待于进一步研究；对于复 杂过渡流型缺乏效果明显的判别方法。 6 、多尺度分析方法 近几年一些研究者开始从多尺度的角度出发，结合非线性理论来研究复杂的多相流 系统。吴浩江等【2 9 j 对管段差压信号利用小波分解一重构方法进行去噪处理后，应用分形 理论提取特定流型对应的特征向量，用于训练人工神经网络，实现了对油水水平管流流 型的自动识别；赵贵兵等1 3 u j 对流化床不同测量位置的压力波动信号进行多尺度分解，并 分别对分解的信号进行r s 分析。研究发现，分解的信号可由多尺度方法得到较好的理 解。 7 、l b m 模拟研究方法 l a t t i c eb o l t z m a n n ( l b m ) 方法是2 0 世纪8 0 年代中期发展起来的一种流场模拟方法， 其理论基础是分子运动论和统计力学。这种方法从模拟流体运动的离散模型出发，在满 足质量、动量和能量守恒的条件下对大量粒子的运动进行统计计算，得到宏观运动参数。 当传统流体力学计算方法不能克服它们存在的不足时，l b m 方法的诞生就受到了人们 高度的关注，希望从不同的角度来研究流体运动的本质问题。在短短的十几年里，人们 对这种方法的理论进行了不断的发展和改进，研究它在不同领域中的应用，这些工作取 第一章绪论 得了非常瞩目的成果，因而本文采用了l b m 模拟研究方法和硬件结合的方法。 1 5 课题主要研究内容及章节安排 本课题采用格子b o l t z m a n n 方法模拟了垂直上升管内的三种流型，采用阵列电导探 针测量方法对垂直上升管内油水两相流进行试验，并进行数据处理和比较，实现了油水 两相流流型的直观识别，为油气水三相流流型研究和水平油井流型的研究打下了坚实的 基础。各章的具体内容安排如下： 第一章对目前油水两相流流动结构及其检测技术现状进行了总结，指出了目前研究 存在的问题及本文的主要研究内容。 第二章介绍了格子b o l t z m a n n 的相关理论知识，如格子b o l t z m a n n 的起源、发展、 优点和应用，重点介绍了格子b o l t z m a n n 的模型和原理。 第三章介绍了油水两相流探针测量系统设计。重点阐述了两相流探针测量方法的测 量原理，介绍了电导探针的设计与制作、测量电路以及信号采集和本次油水两相流的动 态试验。 第四章阐述了格子b l o t z m a n n 边界条件、模型和格子大小等关键问题在本课题中的 选取原则，并给出了当地声速的计算方法，开发了油水两相流流型模拟软件，通过软件 模拟获得了水包油、过渡流型和油包水三种流型的典型模拟图。 第五章分析了垂直上升管油水两相流不同流型电导探针信号特点，对动态试验所测 得的工况进行流型分类，应用递归图对探针信号所做的流型分类进行了验证。确定了动 态试验条件垂直上升管油水两相流流型图，并将其与格子b o l t z m a n n 所得的流型进行了 比较和分析。 第六章为结束语。总结了本文提出的油水两相流流型辨识方法和所取得的成果，并 提出了展望。 第二章l b m 基础知识及其在多相流检测中的应用 第二章格子b o l t z m a n n 基础知识及其在多相流检测中的应用 本章介绍了格子b o l t z m a n n 的相关理论知识，如格子b o l t z m a n n 的起源、发展、优 点和应用，重点介绍了格子b o l t z _ a n a n n 的模型和原理。 2 1l b m 的起源和优点 流体是由流体分子构成的，流体的宏观运动是微观分子热运动的平均结果。因此， 如果能够用计算机模拟流体分子的运动，并能准确地反映出分子之间的相互作用，我们 就可以模拟出流体的宏观运动。这就是分子动力学模拟的基本思想。分子动力学模拟是 一种真正的微观模拟方法 3 1 1 。但是，这种方法的计算代价是巨大的。模拟中，系统演化 的步长必须很小，并且在每一个时间步，需要根据作用在每个分子上的外力和前一时刻 的位置来计算其新位置和新速度。对在前一时刻发生碰撞的任何分子，都要对其进行判 别并计算其新轨道。显然，这需要非常大的计算量。所以，分子动力学模拟目前还只能 用于简单的二维流动模拟p 引，对三维复杂流动进行模拟几乎是不可能的。1 9 8 6 年，法 国的f r i s c h 等人1 3 3 j 在一篇名为“l a t t i c e g a sa u t o m a t af o rt h en a v i e r - s t o k e se q u a t i o n ”的文 章中提出了格子气自动机( l g a ) 的概念。可以说，它就是格子b o l t z m a n n 方法( l b m ) 的 起源。l g a 是一种简化的分子动力学模型。在这种模型中，模拟的对象不是数目庞大 的流体分子个体，而是数目大大减少的流体粒子，即微观充分大、宏观充分小的流体分 子微团。这样，就大大降低了对计算机存储量的要求。可以说，l g a 方法是一种离散 的粒子动力学方法，这是因为它采用的计算网格、时间和速度空间都是离散的。由于 l g a 方法是元胞自动机( c a ) 在流体力学中的具体应用，因此它经常被称为格子气元胞 自动机。 一般情况下，一个格子气自动机包含一个均匀网格和网格节点上的粒子( 即分子微 团) 。用x 表示网格节点，每个节点上有m 个粒子运动方向，可以用一个布尔变量 刀，( x ，) ( i = o ，1 ，2 ，3 ，m ) 表示时刻t 在节点x 上沿某个速度方向是否有一个粒子在运动， 即当有一个粒子运动时，1 1 , = l ；否则，= 0 。l g a 的演化可用如下的方程表示： ，i ( x + 乞，t + 1 ) = ，i ( x ，f ) + f 2 ，( 刀( x ，) )( 2 1 ) 其中e ，是当地粒子速度，q ( n ( x ，f ) ) 称为碰撞算子，它确定粒子接触后如何碰撞。设计 l g a 的碰撞规则时要求碰撞算子满足质量和动量守恒。 q f ( 刀( z ，f ) ) = 0x e , n ，( 刀( x ，f ) ) = 0 ( 2 2 ) 在格子气自动机中，粒子在每个时间步的运动由两个子步构成：( a ) 碰撞：当多个粒子到 达同一个网格节点时，它们按碰撞规则相互作用并改变各自的速度方向；( b ) 流动：每个 粒子按其速度方向移动到最近的网格节点。 第二章l b m 基础知识及其在多相流检测中的应用 第一个l g a 模型是由h a r d y 、d ep a z z i s 和p o m e a u 于2 0 世纪7 0 年代提出的，并根 据三位作者的名字命名为h p p 模型。该模型采用二维正方形格子，每个格点上有流体 粒子驻留。这些粒子具有相同的质量m ，并可沿四个网格线方向e i = ( 1 ，o ) ，e 2 = ( o ，1 ) ， e ，= ( 一1 ，0 ) ，e 4 = ( 0 ，一1 ) 运动。在每个格点上要求最多只有一个粒子沿某个网格线运动， 并在一个时间步内沿网格线运动到相邻的格点。在碰撞和迁移过程中，h p p 模型遵循一 定的物理规律，即保持质量和动量守恒。这正是流体宏观运动方程描述的规律，因此 h p p 模型能够反映出流体运动的一些基本特征。但是，由于正方形格子缺乏足够的对称 性，h p p 模型对应的宏观方程不能反映正确的非线性和耗散效应，即表现为不能导出 n s 方程。 第二个l g a 模型是由f r i s hh a s s l a c h e r 和p o m e a u 3 4 l 于1 9 8 6 年提出的，命名为f h p ， 该模型采用六方形格子。f h p 比h p p 模型有更高的对称性。可以证明，在一定条件下 从f h p 模型可以得到宏观尺度上的不可压n a v i e r - s t o k e s 方程。对称性约束的发现开创 了格子方法研究的一个新局面。 h p p 模型的格子由正六边形构成，每个格点处有六个粒子运动方向： p ，= ( c o s 0 ，) ，p = ( f 一1 ) 万，u = l ，6 ) ，所以每个格点的状态可以由一个位的，sin0, 36 布尔变量表示：刀= ( ，z l ，n 2 ，) ，= 0 或l 表示沿方向e f 运动的粒子数。l g a 模型的 碰撞算子的复杂程度与粒子速度方向的个数有关，一般的m 个粒子速度方向的碰撞算 子的复杂程度为2 肘。图2 1 给出了一些碰撞过程： i n p u ts t a t e o u ts t a t e i n p u ts t a t e o u ts t a t e v 乏 4 * 。v 乏 4 * 4 7 11 5 6 4 7 1 5 6 4 * ， 2 4 * ， lv 一 47 11 5 6 4 * 。v 2 4 * 4 7 1 气 6 图2 1f h p 模型中的碰撞规则 第二章l b m 基础知识及其在多相流检测中的应用 由于采用布尔型变量，l g a 模型不会遇到数值不稳定问题；同时，所有粒子的碰 撞和流动同时发生，且粒子之间的想互作用是局部的，因而l g a 具有本质的并行性， 非常适于在并行计算机上实现；在固壁边界上，流体粒子反弹回流场内部，因而l g a 在处理复杂边晃时非常方便。这些都是l g a 模型不同于传统数值方法的优点。但是， 由于在运行过程中会引入随机噪声，模拟结果往往含有统计噪声。同时，从l g a 导出 的宏观方程与标准的流体力学方程存在一定的差距，如不满足伽利略不变性，并且压力 与速度相关等。虽然通过重新标度可以得到标准的宏观方程，但这种方法只能用于单相 流系统，对于多相流系统不再适用。从计算的角度看，l g a 的碰撞算子具指数复杂性 ( 2 朋) ，需要大的计算量和存储量，这些缺点极大地限制了l g a 方法的应用。 格子b o l t z m a n n 方法就是为了克服l g a 的这些不足而发展起来的。在格子 b o l t z m a n n 方法中，用粒子分布函数的演化代替l g a 中粒子的演化，m c n a n a r a 和 z a n e t t i 3 5 】在1 9 8 8 年首次提出在l g a 中直接使用布尔变量哆的统计平均量厂( 称为单粒 子分布函数) 代替臻进得演化，并且忽略粒子的单个运动特性和粒子间的相关性，最终， 我们就能够得到格子b o l t z m a n n 方程( l b e ) ，即 z ( x + 乞，f + 1 ) = z ( x ，) + q ，( 厂( x ，f ) ) ( f = 0 ，1 ，朋) ( 2 3 ) 在这种模型中，由于f 已经是统计平均后的量，因此模拟结果不再含有统计噪声。 q ，( 厂( x ，) ) 仍是碰撞算子，它表示发生碰撞后，的变化率。q ，( 厂( x ，f ) ) 仅仅与当地分布 函数有关，此外，在l b m 中压力和密度存在这样的关系：p o o p 。但是，这一模型仍然 采用l g a 的碰撞方式，即用q ，( 厂) 代替q ，( 胛) ，所以碰撞算子仍然具有指数复杂性。 h i g u e r e 和j i n e n e z 在1 9 8 9 年对上述模型作了进一步的简化，提出了线性化碰撞算 子模型。该模型引入平衡态分布函数f ) ，对碰撞算子作线性化处理，即 q ，( 厂( x ，嘞= k ，( ，( l f ) 一凹( 墨f ) ) ，其中k ，= a q ，够称为碰撞矩阵。显然，该碰撞算 子的复杂性由2 m 降为m 2 。从数值计算角度看，这种线性化处理简化了计算，并使存储 量大大下降。但是，这种方法的数值稳定性较差。 h i g u e r e 等人在19 8 9 年进一步提出了强化算子格子b o l t z m a n n 模型。在该模型中， 平衡态分布函数的碰撞算子的选择不再依赖原有的l g a 模型，而是根据所要描述的宏 观方程确定。模型的碰撞矩阵是一个对称循环矩阵，且其元素仅仅与离散速度的夹角有 关。 此后，几个不同的研究小组分别独立提出了一种更简单的模型，即单松弛模型。在 该模型中，碰撞过程用趋于某一平衡态的松弛过程代替，矩阵由一个称为松弛时间的参 数确定，即k ，= 一( 1 r ) 。这与分子动理论中b o l t z m a n n 方程的 b g k ( b h a t n a g a r - g r o s s k r o o k ) 算子类似，故这类模型又称为格子b g k ( l b g k ) 模型。即： 1 z ( x + q ，+ 1 ) = z ( x ，f ) 一二 z ( x ，) 一z 叼( x ，f ) 】 ( 2 - 4 ) l b g k 模型极大地简化了模型的计算，并且在一定条件下可以从该模型导出正确的 n a v i e r - s t o k e s 方程，有效地克服了l g a 方法的不足。这类模型的提出使格子方法的研 究达到一个新的水平。目前，此类模型是格子b o l t z m a n n 方法中最主要也是应用最广泛 的模型。 第二章l b m 基础知识及其在多相流检测中的应用 与其它数值计算方法相比，格子b o l t z m a n n 方法具有这样的特点： 1 、l b m 中的流动算子( 或移动过程) 在相空间( 或速度空间) 上是线性的；这个 特性是从分子动力学中继承而来的。在采用宏观量描述的数值计算方法中，对流项是非 线性的。通过多尺度展开技术，可以把结合了碰撞算子的线性流动项转化为宏观的非线 性对流项。在不可压条件下，l b e 能转化成不可压n s 方程。 2 、l b m 中的压力项是采用状态方程计算出来的；相反的，在对不可压n s 方程的 直接数值模拟过程中，压力满足一个p o i s s o n 方程。其中，速度应变是它的源项。通常 在解这样的压力方程时，会遇到一定的困难。这就需要特殊的处理，例如松弛迭代。 3 、l b m 中使用的粒子速度方向个数在相空间上是最少的；在包含 m a x w e l l b o l t z m a n n 平衡分布的传统动力学理论中，相空间是一个彻底的函数空间。它 的平均过程包含整个速度相空间中的所有信息。在l b m 中，由于仅使用了一两个或很 少的几个速度方向，使得微观分布函数和宏观量之间的转换非常简单，因此它的计算算 法也是相当容易的。 4 、l b m 中存在着最适于并行处理的局部相互关系模型； 5 、边界条件容易设定。 2 2l b m 方法的原理及基本模型 l b m 方法是由格子气自动机( l g a ) 发展而来的，因此与l g a 相拟，l b m 方法也可 以看作是一种流体运动的微观模型；另一方面，l b m 可以看作是连续b o l t z m a n n 方程 的一个特殊的离散格式，因此l b m 方法又可以看作是一个介观的动理论格式。下面我 们简单介绍一下格子l b m 方法的原理和其基本模型。 2 2 1l b g k 模型的平衡态分布函数的一般形式及其宏观方程 2 2 1 1 宏观方程 我们从含b g k 模型的连续b o l t z m a n n 方程开始 妥协耵= 一! ( 厂一g ) ( 2 - 5 ) 魂 。 f ” 这里，f 兰厂( z ，c ，t ) 是单粒子分布函数；c 是连续速度空间；f 是松弛时间，g 为局部 m a x w e l l b o l t z m a n n 平衡态分布函数，定义为 咖川= 茄旁e 卅譬， ， 这里，d 是速度空间c 的维数，p 、“和p = k b t m 分别是质量密度、宏观速度及单 位质量的标准化温度，k 、丁和m 分别是b o l z t m n a n 常量、温度及分子质量。质量密 - l o 第二章l b m 基础知识及其在多相流检测中的应用 一一一一一 度p 、速度比和标准化温度p ( 或内能密度) 是厂或g 的流体运动积分 p = 、皿c = 、g d c q _ 7 a ) p “= 肛= f c g d c 争= 硅( c 训2 肚= 硅( ) 2 础 在平衡状态下，f = g ，连续b o l t z m a n n 方程( 2 - 5 ) 变为 丝= o d t ( 2 - 7 b ) ( 2 - 7 c ) ( 2 - 8 ) 旦：旦十c v ( 2 - 9 ) 一= 一十f j d t - o t 通过求解方程( 2 8 ) 的流体运动积分，可以很容易地导出e u l e r 方程 詈柙和炉o ( 2 _ 1 0 a ) _ a ( p u ) + v ( p u u + p o ) = 0 ( 2 - l o b ) a t 丝掣+ v ( 岛+ & + 尸) 甜= o ( 2 - 1 0 c ) 西 、9“ 这里，岛= d p o 2 、s k = p u 2 2 及尸分别是内能密度、平动能密度及压力，“2 = “。甜。 动量方程( 2 1 0 b ) 也可以写为 p 詈w 舻坷尸 ( 2 - 1 1 ) p = p o ( 2 1 2 ) 是理想气体的状态方程。 通过使用c h a p m a n e n s k o g 方法，厂的解可以展开为 f = f o + 厂1 + ( 2 1 3 ) 这里，方程( 2 1 3 ) 中的厂1 要比厂o 小得多。由于o = g ，把方程( 2 1 3 ) 代入方( 2 5 ) f ( 1 ) = - r d j g l l 厂矿- g - f