（流体力学专业论文）球形粒子沉降运动的DLMFD模拟.pdf

独创性声明 本人声明所里交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果，也不包含为获得澎姿盘蓥或其他教育襁 构的学位或证书面使用过的材料。与我一同工作的弼志对本研究所做的任何贡献 均已在沧文中作了明确的说明并表示谢意。 学位论文作者签名：弦好 签字目期：拶g 年占月9f t学位论文作者签名：纭新签字目期：拶占年6 月9 学位论文版权使用授权书 本学位论文作者完全了解澎婆盘鲎有关保留、使用学位论文的规定， 有权保籍并向国家宥关部门或机构送交论文的复印件和磁盘，允许论文被查阅和 借阅。本人授权塑望盘茎可以将学位论文的全部绒部分内容编入有关数据库 进行检索，可以采用影印、缩印或扫描等复制手段傺存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：孑长好导师签名： 讳隧喇 签字日期： 护莎年f 月7 臼 签字日期：? 捌萨孑月夕岛 学位论文作者毕监届去向：成都 工作单位：中国核动力研究设诗院 通讯地址： 电话： 出b 编 浙扛大学硕士学位论文张妍 2 0 0 6 5 摘要 本文在d u l f d 方法的基础上，改变原有的碰撞模型，分别采用硬球模型和 软球解析模型来处理球体之间以及球体和壁面之问的作用，数值模拟了球形粒子 在三维方槽中的沉降。通过三种不同碰撞模型的对比，研究了水动力和短程力在 粒子浓度不同时的主导作用。并采用软球解析模型研究了雷诺数、初始位置等几 个重要参数对球形粒子沉降的影响。 在对球形粒子沉降的研究中，首先采用硬球碰撞模型和软球解析模型进行模 拟，得到的结果与实验非常吻合，从而验证了两种模型的可行性。当沉降粒子群 浓度非常低如只有单个或两个圆球参与沉降时，对三种模型模拟得到的轨迹，速 度进行分析得出结论：无论采用何种模型来处理碰撞的作用，得到的结果是非常 接近的，也就是说在圆球并不密集的沉降中，对球体沉降影响较大的是其中的水 动力，而不是碰撞时的短程力，完全可以从计算简便快速的角度选择其中一种模 型。而对2 5 个球形粒子的沉降模拟结果进行分析，可以发现，由予处理圆球之 间以及圆球和壁面之间作用的方式不同，三种模型得出的沉降轨迹、平均沉降速 度等不再相似，说明粒子群浓度较大时，起主导作用的是短程力而不再是水动力。 对于三种不同的模型，对碰撞作用的处理方式有很大的不同。原有的d 瑚f d 方法中采用的是简单排斥力模型，粒子不能真实接触发生碰撞，其中涉及的排斥 力不仅和球心位置有关，还和一些碰撞参数有关，而这些碰撞参数并非定值，选 择的数值不同，结果也有变化；软球解析模型在低雷诺数时有明显的物理意义， 粒子碰撞时能真实接触甚至有些许重叠，碰撞时的作用力和间距、粒子速度以及 粒子本身的物理性质有关，不受其它数值不定的参数影响，相对排斥力模型有一 定的优势；硬球碰撞模型从能量的角度出发，假设粒子碰撞时能量守恒，碰撞前 后直接改变粒子的速度，方法简单、计算简便。采用软球解析模型模拟了不同雷 诺数下单个、两个以及2 5 个圆球在三维方槽中的沉降。在单个圆球沉降轨迹的 x y 视图上，低雷诺数呈一条与壁面平行的线段，而高雷诺数情况下则在平行的 两个壁面之间来回振荡或者沿四个壁面依次碰撞。在两个圆球的沉降中，初始位 置的改变对低雷诺影响较大，如初始两个圆球肩并肩释放将一直维持这种状况向 下运动，但在高雷诺数时则不是如此，雷诺数的不同还会改变两个圆球的球心距。 2 5 个圆球沉降中，初始均匀放置，和二维时多个圆形粒子沉降的现象相似，能 3 苎塑翌苎型竺丝苎 垡塑! ! 塑：! 观察到粒子整体由对称结构一一v 结构一一混乱，观察到了球形粒子的粒子群现 象 关键字：d l m f d 方法；简单排斥力模型；硬球碰撞模型；软球解析模型； 球形粒子的沉降；雷诺数；初始位置； 4 浙江大学硬士学位论文张妍2 d 0 6 5 a b s t r a c t t h ed l m f dm e t h o di sa p p l i e dt ot h es e d i m e n t a t i o no fs p h e r e si nt l l 阳e - d i m e n s i o n a lv e r t i c a lt u b e t h eb i n a r yh a r dm o d e la n ds o f t - s p h e r em o d e la r e i n t r o d u c e di nt h ec o l l i s i o n 咖a t e g yt od e a lw i t ht h ei n t e r a c t i o n si nc o n t a c t t w on e w m o d e l sa 北c o m p a r e dw i t ht h es i m p l er e p u l s i v ef o r c em o d e l a n dt h ed o m i n a n tr o l eo f l o n g - r a n g eh y d r o d y n a m i cf o r c ea n ds h o r t - r a n g ei n t e r a c t i o n s a r ee s t a b l i s h e di n d i f f e r e n tc a s e s e f f e c to fs o m ep a r a m e t e r sl i k er e y n o l d sn u m b e r , i n i t i a lp o s i t i o n sa r c s t u d i e du s i n gs o f t - s p h e r em o d e lu n d e rd i f f e r e n tc o n d i t i o n s f i r s t , e x p e r i m e n t a lr e s u l tv a l i d a t e dt h eh a r d - s p h e r em o d e la n ds o f t - s p h e r em o d e l t h r o ec o l l i s i o nm o d e l sa r ci i i t r o d u c e dt os i m u l a t et h em o t i o no f s p h e r e si ns u s p e n s i o n , a n dt h ea n a l y s e so fw a k e sa n dv e l o c i t i e ss h o wt h a t ：i nan o n - c o n c e n t r a t e ds u s p e n s i o n c o n t a i n i n go n l yo n es p h e r eo rt w os p h e r e s ，t h em o t i o ni si n d e p e n d e n to ft h ec o l l i s i o n m o d e lc h o s e ni x 娥l u s et h eb e h a v i o ro fs p h e r e sh e r ei s m a i n l yc o n t r o l l e db yt h e l o n g - r a n g eh y d r o d y n a m i cf o r c e s ，a n dt h ee f f e c to f t h es h o r t - r a n g ei n t e r a c t i o n st a k i n g p l a c ea t n e a rc o n t a c ti s i n s i g n i f i c a n t b yc o n t r a s t , r e s u l t sa r cd i f f e r e n tf o rt h e c o n c e n t r a t e ds u s p e n s i o n ：t h es e , d i m e n t a f i o no f2 5s p h e r e si nat h r e e - d i m e n s i o n a l v e r t i c a lt u b ea g es i m u l a t e du s i n gt h r e ed i f f e r e n tc o l l i s i o nm o d e l s ，a n dt h ev e l o c i t i e s a n dw a k e sa r en o ts i m i l a ra n ym o r e t h er e s u l td e m o n s t r a t e st h a ti nac o n c e n t r a t e d s u s p e n s i o n , t h es h o r t - r a n g ei n t e r a c t i o n sa f f e c t t h em o t i o ns i g n i f i c a n d y t h r e ec o l l i s i o nm o d e l sa r eb a s e do nd i f f e r e n tp r i n c i p l e s t h ef u n c t i o no fs i m p l e r e p u l s i v ef o r c eu s e di nt h eo r i g i n a ld l m f dm e t h o di sn o to n l yr e l a t e dt ot h e d i s t a n c eo ft w os p h e r e sb u ta l s oa f f e c t e db yt h ec o l l i s i o np a r a m e t e r sw h i c ha 托n o t p h y s i c a l f o rt h es o f t - s p h e r em o d e l , 缸t h el o w e rr e y n o l d sn u m b e r , t h ef u n c t i o no f s h o r t - r a n g ef o r c eh a sp h y s i c a lm e a n i n g ，a n dt h ef o r c ei so n l yr e l a t e dt ot h es e p a r a t i o n a n dt h e p a n i c l ev e l o c i t i e s c o m p a r e dt ot h er e p u l s i v ef o r c em o d e l ，r o t a t i o n i s c o n s i d e r e dw h e nt h ec o l l i s i o nh a p p e n sa n dt h ep a r t i c l e st o u c ho re v e no v e r l a p f o rt h e h a r d - s p h e r em o d e l ，t h em o m e n t u me x c h a n g eb c t w e e nt w oc o l l i d i n gp a r t i c l e si s c o n s i d e r e di n s t e a do ft h ef o r c e s p o s t - c o l l i s i o nv e l o c i t i e so ft h ep a r t i c l e sa r e d e t e r m i n e df r o mt h es y s t e mm o m e n t u ma n de n e r g yc o n s e r v a t i o no rl o s s ，w h i c hb r i n g s 5 浙让大学硕士学位论文 张妍 2 0 0 6 5 n u m e r i c a lc o n v e n i e n c e i nt h i st h e s i s , 啪s u p p o s ec o l l i s i o n sa r ea l le l a s t i c i nt h ef o r t h p a r t , v a r i o u sc a s e sa r es i m u l a t e dw i t ht h es o f t - s p h e r em o d e l ，a n dt h er e s u l t sf o r d i f f e r e n tr e y n o l d sn u m b e r s ，d i f f e r e n tn u m b e ro fs p h e r e sa n di n i t i a lp o s i t i o n sa r e g i v e n c o n c l u s i o n sa r i e ：a tl o w e rr e y n o l d sn u m b e r s ，t h em o r es p h e r e se x i s t ，t h eh i g h e r a v e r a g es e t o i n gv e l o c i t i e sa r e f o l t h es e d i m e n t a t i o no f o n es p h e r e ，t r a j e c t o r i e so f t h e s p h e r ei nx y = p l a n ev i e wa td i f f e r e n tr e y n o l d sn u m b e r sf i r et y p i c a l ：t h es p h e r em o v e s a l o n gt h el i n ep a r a l l e l e dt ot h ea x i sa tl o w e rr e y n o l d sn u m b e r s a n do s c i l l a t e sa r o u n d t h et w op a r a l l e lw a l l so rt o u c ha l lw a l l sc l o c k w i s ea th i g h e rr e y n o l d sn u m b e r s f o r t h es e d i m e n t a t i o no f t w os p h e r e s ，t h ee f f e c to fi n i t i a lp o s i t i o n si si m p o r t a n tf o rl o w e r r e y n o l d sn u m b e r s ：w h e nt w os p h e r e sa r cr e l e a s e ds i d eb ys i d e ，t h e yk e e pt h es t a t ea l l o v e rt h es e d i m e n t a t i o np r o c e s sa n dn oc o l l i s i o nh a p p e n s , w h i l er e l e a s i n ge n d - t o - e n d v e r t i c a l l y , ap h e n o m e n o nc a l l e dd k tc a nb eo b s e r v e d b u ta th i g h e rr e y n o l d s n u m b e r s , t h ei n i t i a lp o s i t i o n sa f f e c tn e g l i g i b l y i nt h es e d i m e n t a t i o no f 2 5s p h e r e s ，t h e v a r i e t yr u l eo ft h ec o n f o r m a t i o no ft h es p h e r e se n s e m b l ei ss i m i l a rt ot h ec i r c u l a r p a r t i c l e si nt w o - d i m e n s i o n a lc h a n n e l , ar a n g eo ft y p e so fs t a t e sl i k es y m m e t r y , “v ” $ t l u c t u r e $ a n do u t - o f - o r d e ra tt h el a s tc a nb eo b s e r v e d ，a n dt h et y p i c a lp h e n o m e n ao f c l u s t e r i sa l s oi n t r o d u c e di nd e t a i l k e y w o r d s ：d l m f dm e t h o d ；h a r d - s p h e r em o d e l ；s o f t s p h e r em o d e l ；t h e s e d i m e n t a t i o no f s p h e r e s ；r e y n o l d sn u m b e r s ；t h ei n i t i a lp o s i t i o n s ； 6 浙江大学硕士学位论文 第一章概述 第一章绪，论 1 1 论文的背景及理论概述 1 1 1 粒子的沉降 粒子悬浮流通常指固粒包含在流体中的流动，它在自然界和工业中的应用非 常广泛，如江河中的泥沙沉积流化床中的气固流动等，还被普遍应用到粉尘洁 净、造纸过程以及食品加工等实际应用领域。而固体粒子在重力作用下的沉降运 动作为粒子悬浮流的一种情况，应用也相当普遍，因为尺度和形状迥异的粒子在 雷诺数不同的情况下具有不同的沉降运动特性，这种特性上的差异可以用来进行 粒子和流体之间或者不同粒子之间的分离，例如水煤浆的脱水，废水的净化处理 等等。本文重点研究球形粒子的沉降，很多情况能直接推广到对异形粒子沉降的 研究。 对球形粒子沉降的流场研究由来已久，在日常生活中，往往出现的是大量粒 子共同沉降的过程，这其中粒子与粒子之间以及粒子和外壁之间的碰撞非常频 繁，有结论表明，两相流中圆球粒子沉降的平均速度会因为粒子间的相互作用以 及流体对粒子的水动力作用而减慢。一般说来，固粒在流体中沉降主要受到流体 的浮力以及粘滞力的作用，而每个固粒i 的平动和转动都要受到周围其它粒子的 影响，主要表现在两个方面，一为l o n g - r a n g eh y d r o d y n a m i ci n t e r a c t i o n s ，称之为 水动力作用，是其它粒子通过流体施加给粒子i 的，它会随着两个固粒之间距离 的增加而很快衰减，也称为长距影响；对应长距的为短距的作用，称为s h o r t - r a n g e i n t e r a c t i o n s ，通常包含润滑力( 1 u b r i c a t i o nf o r c e ，又称短程水动力) 以及为防止润 滑力作用下会导致粒子重叠而产生的短程排斥力( 也称接触力) 。粒子沉降过程 中的这些相互作用都非常的复杂，在实际的直接数值模拟中，往往通过解析的式 子对它们进行简化或者对应具体问题忽略若干作用进行模拟。 在数值模拟粒子的沉降时，粒子间以及粒子和外壁间的碰撞是一个难点，现 有的理论方法都很难较形象地反映出这一瞬间，因为真实的碰撞与雷诺数、材料 性质以及表面粗糙度等都有关。c c r o w e 1 j 等( 1 9 9 8 ) 通过研究得出结论：一般说来， 7 浙江大学硕士学位论文 第一章概述 粒子悬浮流中的碰撞模型可分为两类：硬球碰撞模型( b i n a r yh a r d - s p h e r em o d e l ) 和软球模型( s o f t - s p h e r em o d e l ) 在硬球碰撞模型中，认为两个粒子接触的瞬间 就交换动量，瞬间碰撞后的速度取决于系统的动量( 不) 守恒和能量( 不) 守恒， 也就对应了通常所说的( q t e ) 完全弹性碰撞两种情形。而软球碰撞模型，则认为 发生碰撞的两个粒子的速度由牛顿运动方程来决定，其中的作用力是可能与两个 粒子间的距离甚至是粒子碰撞前的速度相关的函数。在粒子并不密集的 ( n o n - c o n c e n t r a t e d ) 悬浮流中，粒子的行为大多受到上面提到的水动力作用，而 受短程作用不大，因此这种情况下，只需从计算简便的角度任意选择一种模型来 进行数值模拟。粒子密集的情况，软球碰撞模型应用较多。此外，在粒子雷诺数 很低( 粒子惯性也很小) 时，如球形粒子的表面光滑，短程水动力( 1 u b r i c a t i o n f o r c e ) 理论上可以阻止粒子问的碰撞；如果表面并不光滑，则粒子间仍可接触 从而产生排斥力在粒子惯性较大时，会出现较强的粒子碰撞，在数值模拟中可 采用硬球非碰撞或软球碰撞模型( 粒子浓度很大时) 。凡此种种，都为了尽可能 使数值模拟的碰撞模型反映真实的物理碰撞过程。 i 1 2d l n f d 方法的发展 根据作用的水动力方式的不同，将近二十年内涌现的多种水动力模拟方法分 为两类：直接和非直接数值模拟方法，分别简称为d n s 和n o n d n s 方法。其中 在1 9 9 2 年h u 埘，u n v e r d 3 1 4 1 首次提出多相流的直接数值模拟方法。直接数值模拟 和非直接数值模拟方法的主要区别在于：前者是直接利用n - s 方程来计算水动 力，而后者则不是。常见的n o n - d n s 方法有p a r t i c l e b a s e dm e t h o d s 包括m o l e c u l a r d y n a m i c s 方法似d ) 和d i s s i p a t i v ep a r t i c l ed y n a m i c s 方法( d p d ) 等，还有其它如 f t s 、s d 方法等等。而d n s 方法【坷【6 】通常涉及到有限元方法、边界元方法、有 限差分法、控制体积法和格子玻尔兹曼方法等来计算流场，一般依据计算流场时 更新网格的方式分成两类：基于贴体网格的方法和基于非贴体网格的方法，前者 在流体占据的物理空问上求解n s 方程，而后者在求解时，流场的计算区域不 仅包含流体占据的空间还包括了固粒所占的空间，主要包括格子一玻尔兹曼方法 ( 1 a t t i c e b o l t z m a n nm e t h o d ) 1 7 】和虚拟区域方法( f i c t i t i o u sd o m a i nm e t h o d ) ，后来 虚拟区域法又渐渐发展成虚拟边界法( b o u n d a r ym e t h o d ) 和基于分布式拉格朗 浙江大学硕士学位论文 第一章概述 日乘子的虚拟区域方法( d i s t r i b u t e dl a g r a n g em u l t i p l i e rb a s e df i c t i t i o u sd o m a i n m e t h o d ，以下简称d l m f d ) 等。 本文用来求解粒子沉降问题的主要方法则是上述所说的d l m f d 方法，它 是在1 9 9 4 年由6 1 0 w i n s k i 等人嘲从虚拟区域法演化得到的。虚拟区域法中，计算 边界包括原始的边界和通过施加边界约束条件得到的虚拟边界，使得计算区域更 加简单明了，而o l o w i n s l d 等人嘲f 9 1 1 1 0 l ( 1 9 9 4 ) 在d i r i c h l e t 问题上通过对边界施 加拉格朗日乘子进行约束，诠释了一种新的虚拟区域法即d l m f d 方法。起初， d l m p d 法用于微分方程的求解，g l o w i m l d 等人将它应用到非线性的时间相关 河题，如管道中的粘塑性介质的流动、g i n z b u r g - l a n d a u 方程和n s 方程的求解 等等。后来，b e r t r a n d 等人( 1 9 9 7 ) 和t a n g u y 等人( 1 9 9 6 ) 用它来计算三维牛 顿流体的s t o k e s 流以及混合容器中粘一塑性流体的流动问题。1 9 9 7 年， l o w i n s k i 等人【1 1 1 又通过在边界上施加拉格朗日乘子实现无滑移边界条件，从而成功模拟了 绕某一有确定运动的粒子的非定常粘性不可压流动。1 9 9 8 年【1 2 1 ，他们再次选择 了合适的拉格朗日乘子对边界施加约束，实现了粒子内部的刚化，对计算区域模 拟的结果同粒子边界上仅有无滑移边界条件约束的情形非常吻合。d l m f d 方法 用于粒子悬浮流的研究始于1 9 9 7 年，g l o w i m k i 等人 1 3 1 在1 9 9 9 年又引入了粒子 间以及粒子和外壁问碰撞时的简单短程排斥力模型，并完整地描述了该方法，标 志着用于直接数值模拟粒子悬浮流的d l m 方法的出现。 d l m f d 方法的基本思想可以分两个方面来理解，一方面计算区域不仅包含 流体运动区域还包括被粒子占据的那部分区域，并且后者的密度认为和流体密度 相当，这样一来就将计算区域变得规则，称为f i c t i o u sd o m a i n ；另一方面，分 布式拉格朗日乘子则作用在固粒表面将粒子内被流体占据的区域，从而满足刚体 运动的约束条件。d l m f d 方法有如下几个要点：( 1 ) 公式基于守恒型的方程， 便于处理非牛顿流体；( 2 ) 对于异形粒子，边界条件相对复杂，计算也可在简单 区域内进行，而边界条件只需调节拉格朗日乘予来实现制约；( 3 ) 不需要通过显 示计算粒子上的力来确定粒子的运动，因此计算可以更加稳定。 s m g h 等人0 4 1 ( 2 0 0 3 ) 将该方法进行扩展并应用到粘弹性流动中。同年，发 展了另一个版本的d l m 方法( p a n t a n k a r 等人，2 0 0 0 ) ，解决了o l o w i n s k i 版本 中不适用予中性悬浮粒子的情况，然而由于其中的拉格朗目乘子必须定义在科 9 浙江大学硕士学位论文 第一章概述 空间上，执行起来非常麻烦最终没能继续发展。在随后的几年里，p a n ， j o s e p h & g l o w i n s k i ( 2 0 0 1 ) 应用d l m f d 方法，模拟了二维空腔中的6 4 0 0 个圆形 粒子在泰勒不稳定性下的沉降，p a n 等人【1 6 1 ( 2 0 0 2 ) 还研究了1 0 2 4 个球形粒子 的流化问题，这些模拟结果都和实验非常吻合。同一时期，d l m f d 方法还用于 对异形粒子沉降的研究，如g l o w i n s k i 等人( 2 0 0 1 ) 、j u a r e z 等人【1 7 1 ( 2 0 0 2 ) 模拟 t - - 极形粒子的沉降，p a n 等人( 2 0 0 2 ) 1 8 1 模拟了牛顿流体中的椭球沉降等等。另 外本文的指导老师邵雪明教授也应用d l m f d 方法对不同霄诺数下单个圆形粒 子在二维垂直方槽中的沉降特性进行了研究和相关的实验验证【嘲；同时还对两 个相同刚和不同大小圆形粒子之间的相互作用【2 1 1 ，6 0 0 个粒子的沉降过程l ，方 形粒子的沉降过程 2 3 1 ，以及粒子在o l d r o y d - b 流体中的沉降 2 4 1 都做了研究。2 0 0 2 年，y u 等人【2 咒用“q i - p o ”有限元格式代替g l o w i n s k i 等人采用的一阶精度的算 子分裂格式和采用三角形单元离散的有限元格式，使计算精度大大提高。2 0 0 4 年，y u 等人瞄l 又首次采用基于半交错网格的有限差分法取代了原来的有限元方 法，大大提高了方法的计算效率，使得d l m f d 方法迈入了一个新的里程碑。 1 2 国内外研究现状综述 1 2 1 对球形粒子沉降的研究 球形作为粒子形状的一种典型形式，在过去的实验和数值模拟中很早就有了 细致的研究。c h r i s t o p h e r s o n 和d o w s o n 鲫( 1 9 5 9 ) 研究了圆球在椭圆容器中的 沉降特性，指出球体旋转沉降并偏离中心线。t a n n e r t 船l 进一步发现在非牛顿流 体中这种现象更加强烈，并指出在沉降中可能存在圆球偏离中心和旋转的临界球 半径。l e a l 【捌研究了在s t o e k e s 流中圆球粒子的沉降问题。g o l d m a n 等人p o l 论述 了圆球和壁面的吸引作用，同时也研究了圆球沿着倾斜的壁面滚动并指出它们不 是真正的接触，丽是滑动，当圆球被迫靠近壁面时将产生不规则旋转b u n g a y 和 b r e n n e r 【3 1 1 给出了球沿着垂直管壁改变方向的情况，这时两者的距离趋向于零。 h u m p h r e y 和m u r a t a 指出当管倾斜角逐渐增大或者球接近壁面对，旋转的球从正 ( 反向下坡旋转) 逐渐变成负( 下坡旋转) 。g l o w i n s k i 等人l j 3 】( 1 9 9 9 ) 应用d l m f d 方法研究粒子悬浮流，验证了两个球体完整的d k t ( d r a f t i n g ，k i s s i n g ， l o 浙江大学硕士学位论文第一章概述 t u m b l i n g ) 过程。y u 等人陶( 2 0 0 4 ) 应用g l o w i n s k i 等人的方法，模拟了圆形容 器中，圆球在不同的雷诺数下的沉降特性，观察到了高雷诺数情形下经典的发夹 涡现象。 1 2 2 碰撞模型的改进 对于数值模拟过程中出现的球体与球体的碰撞问题，直接数值模拟发展的早 期，人们人为地引入粒子之间的排斥力，一般在粒子表面间距小到一个速度网格 步长时起作用，从而避免碰撞发生。标准的d l m f d 方法中，g l o w i n s k i ，p a n ， h e l s a ，j o s e p h 以及s i n g h ，j o s e p h ，h e s l a ，g l o w i n s k i ，p a n 【1 4 】等还要确保乘 子网格与速度网格不重合，这些必然带来方程的超正定情况。在粒子悬浮流中， 研究粒子之间的碰撞或者说n e a r - c o l l i s i o n s 存在非常多的困难。实际上，连续 介质系统中，表面光滑的粒子不可能在有限的时间内发生碰撞。然而，在某些模 拟中，粒子间距异常的小，加上数值模拟的截断误差，可能导致模拟中发生碰撞 甚至粒子重叠的情况。研究粒子与粒子之间的碰撞时，由于发生环境是两个粒子 间距非常小的情况，为了精确，势必要求网格更加细化，计算精度高，花费也就 相应变大。在应用d l m f d 方法进行数值模拟时，g l o w i n s k i 等人【1 3 】( 1 9 9 9 ) 引 入了短程排斥力的简单碰撞模型，认为粒子间距达到某一指定数值时产生排斥 力，且随着粒子间距的减小而不断增大，直至两粒子分离，在这个模型中，由于 给定数值和网格方面的原因，粒子间始终保持一定距离，无法真正接触。s i n g h 等人 t 4 l ( 2 0 0 3 ) 在前者的基础上改进碰撞模型，允许粒予间任意靠近甚至略微 重叠，他同样引入了上述的排斥力碰撞模型，当粒子重叠对，他的做法是：在运 动方程中将剐体运动限制部分冲突的某一项简单地除去一一这一项可能正是对 应粒子表面上距离结点较远的项。这样做，能避免方程组的超正定情况，但是仅 仅局限在两粒子问重合的距离不会超过一个单元格的情形。但g l o w i n s k i 和 s i n g h 两个模型都没有考虑碰撞对粒子转动的影响。 1 3 本文的研究内容及主要创新点 本文主要是在g l o w i n s k i 等人提出的d l m f d 方法的基础上，研究改进碰撞 模型后模拟一个和多个球形粒子在牛顿流体中的沉降。 创新点主要有： 羔塑苎兰燮竺堡壅 苎二兰 堡垄 ( 1 ) 碰撞模型采用了硬球模型： ( 2 ) 将格子玻尔兹曼方法中介绍的软球解析碰撞模型也引入了d l m f d 中。 浙江大学硬士学位论文 第二章基本理论和方法 第二章基本理论和方法 本文研究球形粒子在方槽中的沉降运动时，采用9 l m f d 方法计算流场，用 硬球模型和软球解析模型分别模拟粒子之间以及粒子和外壁之间碰撞时的相互 作用，以下将分别介绍这些理论方法。 2 1 基于分布式拉格朗e l 乘子的虚拟区域方法简介 2 1 1 组合动量方程 图1 计算区域示意图 图l 所示是本文中圆球在三维方槽中沉降的计算区域示意图，其中只代表球 形粒子，圆球粒子的沉降在图示计算区域为0 的方槽中进行，q 既包括流体占 据的方槽空间，也包括粒子只内部和边界所占据的空间，r 代表区域q 的所有外 边界，为了以后说明的方便，我们定义只特指第i 个球形粒子的内部，而用嬲专 指第i 个粒子的边界部分。研究粒子沉降运动的d l m f d 方法主要由以下的流体 和刚体控制方程组成基本计算公式： 1 3 兰鬯童曼塑主型苎 兰三烹苎查墨丝塑查鲨 ( 1 ) 流体控制方程 ( 2 ) 刚体控制方程 ( 3 ) 动力学方程 见害= 见寥丹盯加q 、雨 v u 2 0 nq 、即) u 。i i r ( f ) o nf = u ，+ ，毛 o na 只( f ) i - - l ，n u 【- o = o nx p ( o ) m d u g ：m g + e d t i”l 业立：互dt 1 u ，l = i j j 4 几- o = 棚 坚：i j 。 d t 堕：， d t x , l 。= k 。 ( 2 1 ) ( 2 2 ) ( 2 3 ) ( 2 4 ) ( 2 5 ) ( 2 6 ) ( 2 7 ) ( 2 8 ) ( 2 9 ) ( 2 1 0 ) ( 2 1 1 ) ( 2 1 2 ) o 见0 = e 柚 ( 2 1 3 ) 其中， 粤：粤+ ( u v ) u ( 2 1 4 ) 毋西、 一 指的是随体导数； 尸( f ) = u p s ( t ) 表示固粒占据的区域；弓- - - - i k ；其它参数见，h 和。分别代表流体的密度， 速度和应力；m ，正，u 。，q ，墨和0 ，分别表示第i 个粒子的质量，惯量， 平动速度，角速度，质心位置和角位移；此外， 1 4 ii-口，li， 浙扛大学硕士学位论文 第二章基本理论和方法 e = l ( f ) m l d s ，霉= l ( ，) 巧x m l d s ( 2 1 5 ) 指的是作用在第i 个粒子上的水动力和相应的力矩，五指由粒子边界指向粒子外 的单元法向量。初始速度值u o 和边界速度值蜥o 满足如下的兼容性条件： v u 。= o nq 、而h ( f ) r i d s = 0 ( 2 1 6 ) 本文中的流体均认为是牛顿流体，应力o = - p i + 2 r d u 】。 ( 2 1 ) 一( 2 1 6 ) 再根据d i r i c h l e t 边界条件可以得到组合动量方程如下： 对于f ) ，o ，寻找，e 露( 固，a a ( d ，u e r 2 和r 满足： l e + ( u d u 】吡一l 胛吡+ 面2l d 【i i 】：d 【v 协+ ( 纠) m 哿一乃护+ 丛铲嘲- v v - - t , 弘( 加_ ( v 吩 所有的v ，v c r 2 和亏r 2 ( 2 。1 7 ) l q v u d x = 0 对所有的q f ( 囝 ( 2 1 8 p ( o = o 对所有的a ( f ) ( 2 1 9 ) 其中，4 指固体粒子和流体的密度比：p 和五分别表示压力和l a g r a n g e 乘子； v ，v ，孝，q 和则是相应的变量。上述的空间分别定义如下： = v f ( q ) 2 l v = i l r ( r ) o n f ( 2 2 0 ) = v ( q ) 2 l v = oo n r ( 2 2 1 ) 露( 锄= g f i l 啦= o ) ( 2 2 2 ) 2 1 2 计算方法 为了求解上节中提到的公式( 2 1 7 2 1 9 ) ，g l o w i n s k i 等人提出算子分裂 格式并最终将方程分解成几个子系统问题进行求解。主要的思想则是将式子中涉 及到的流固系统部分分成相应的流体运动和刚体运动部分， 1 。流动部分：寻找u ”1 彪e 辟0 ，p e 露( 回满足： 1 5 浙江大学硕士学位论文 第二章基本理论和方法 l f i ! n + 1 1 石2 一u n + ( d u 】v 叔一l 胛v 出+ 素l d 【u 】：d 【v 】斑= 。( 2 2 3 ) 所有的v w o 【i 可v n 州总斑- - - 0 对所有的q 工2 ( q ) ( 2 2 4 ) 式子( 2 2 3 2 2 4 ) 代表n - s 方程和连续性方程的弱解形式，可表示如下： 百g u + u v l i = - 唧+ 去v 2 u ( 2 2 5 ) v - u = 0 ( 2 2 6 ) 早期，人们( g l o w s k i n t n 】，1 9 9 7 b ) 采用有限元方法( 如图2 ( a ) 所示) ， 用三角形单元格来离散区域进行求解，后来发现这样做导致粒子速度相对中心线 不是对称的，于是改采用矩形单元格来离散，优点主要体现在这种离散情况下的 对称性相当好。本文中我们则采用有限差分法( 如图2 ( b ) 所示) 近似求解来 取代有限元方法。一般说来，对规则区域有限差分法要比有限元方法效率高的多。 r n 弋溉n 。 i r 乐呵n 八八 ，乐v n 除。v l ( a )( b ) 图2 ( a ) 有限元法三角形单元离散；( b ) 半交错网格。其中表示压力离散点 在采用上图所示的半交错网格离散计算区域后，式( 2 2 3 2 2 4 ) 可以利用 一 常用的投影法进行更进一步的划分，便于求解： ( i ) ( i i ) 1 u * - - r u n + ( u 丽) 4 = 去v 2 u 其中n = u r 伽r ( 2 2 7 ) 一ll，iiil， 一训一叫一训一州叫一州一 盯iilili一 一刮训一州一叫一卅一训一一叫叫一叫一叫叫一叫一 一ili一一ii一 一训一训一训一州一训一训一 一叫一训一训一叫叫州一 叫一州一训州一叫一训一 浙江大学硕士学位论文 第二章基本理论和方法 v 2 p ：宰其中罢：o 。一r f一锄 ( i i i ) i i n + l 石2 - - 一u * = 唧 式子( 2 2 5 ) 若不考虑粘弹性项则可认为是一个对流扩散问题， 方法分成三步得到；4 - - 对角矩阵进行求解，一般是如下形式： uk+13_uk+(uvu广：上2uk+113-一，2产k+vy)，3rf 、。 e 、4 7 。 u k + 2 面1 3 _ _ u 广k l 3 + ( u 帆尸= 去v ；u k + 2 1 3 2 3 ) ， 百u f + l n k + 2 1 3 + 似肌) 删3 = 去r v ；u k + 2 3 v 2 ) 其中。我们设宦一；i i n u ；n k + 1 ( 2 2 8 ) ( 2 2 9 ) 采用a d i ( 2 3 0 ) 式子( 2 2 6 ) 为一典型椭圆方程，诺曼边界条件下采用快速余弦变换( f c t ) 法来求解c s 钔嘲城义函数五：篓细掣及其反函数 乃= 艺k = o 最c o s 业n ，即演变成对不同情况下椭圆方程v 2 = g 的求解，文中 不再详细叙述。 利用( 2 2 5 ) 和( 2 2 6 ) 分别求得h 和p ，代入式( 2 2 7 ) 中可得一，如 此循环计算最终解决了组合动量方程的流体部分。 2 。粒子运动部分： - , 寻找u c j ，u e r 2 ，m r ，九a a 满足： 唐咄+ m e 一乃争v + 亟司一 ( 2 3 。) f - v - x , 亏= ( 五，v 一( v + 号r ) ) ，( ，)其中v e ，v e r 2 ，写e 震 掣，u 一( u + 国r ) r = o 其中e a ( 2 3 2 ) 1 7 塑坚盔兰堕