（理论物理专业论文）晶格玻尔兹曼方法对血液流的初步研究.pdf

晶格玻尔兹曼方法对血液流的初步研究 摘要 自从1 9 9 0 年基于单弛豫时间近似的晶格玻尔兹曼方法( 1 a t t i c eb o l t z m a n n m e t h o d ，简写为l b m ) 提出以来，对l b m 方法及其应用的研究引起了人 们广泛的兴趣，关于l b m 及其应用每年一届的国际会议已连续开了十二 届。目前l b m 已被普遍认为可以作为一种新的流体力学计算方法，而其高 效、精确和鲁棒性( r o b u s t n e s s ) 已得到广泛的证实。该模型已被广泛应用于 研究湍流、两相流、反应扩散系统、颗粒流和悬浮体等等系统。l b m 是一种 分立( d i s c r e t e ) 统计模型，计算过程中信息传递是局域的，因此并行运算程 序关联度很小，各个n o d e 之问需要传递的信息量很小，特别适合在p c 机 群( p c c l u s t e r ) 甚至i n t e r n e t 网络上进行并行计算。用l b m 方法可以在费 用比较小的情况下进行较大规模的数值模拟。血液流系统是有代表性的边界 特别复杂的生物系统。本文中我们用l b m 对该系统的大小动脉内的血液流 进行初步探索。我们的工作主要包括以下几个方面 ( 1 )建立了存在可形变边界的情况下，流体对边界的作用力的精 确计算方法。由于在晶格玻尔兹曼方法中，粒子的微观速度是离散的。而且 为了提高计算效率和减少复杂性，速度集里速度的个数在满足对称性的条件 下尽可能越少越好一般在二维的情况下也就几个速度，三维情况下_ 卜几个 速度。如此一来，当边界是与网格主线斜交的直线( 三维是平面) 或者是曲 线( 三维是曲面) 时，用边界所获得的净动量来计算边界所受到流体的作用 力将出现很大的噪音。虽然流体在边界上的作用力在长直边界上的平均是准 束鳕押p 雳、导师硒薅 钶垒文公布 确的，但是对局部力的计算的不准确有可能导致可形变边界的形变失控。我 们建立了新的基于压力张量积分的方法来计算流体对边界的作用力，可以很 好地解决这个个问题。我们模拟了单个粒子在流体中的沉降，流体对边界的 作用力分别用压力张量积分法和动量交换法来计算。和有限元的计算结果比 较，用压力张量积分法算得的结果符合得相当精确而用动量交换法算得的粒 子在水平方向的速度符合得不太理想。在p o i s e u i l l e 流中，当悬浮粒子的直 径和管道直径之比d s d o 2 时，悬浮粒子会迁移到偏离中心线的位置，这 个现象叫s e g r 6 - s i l b e r b e r g 效应。压力张量积分法能够正确反映该效应而动 量交换法却不能。 ( 2 )建立了用l b m 数值模拟包含有红血球的血液在有狭窄的小 血管中流动的二维模型。为了降低复杂性，参照文献【4 7 ，4 8 的做法，这犁 红细胞是不可形变的圆柱。狭窄由两个对称的半圆柱构成的鼓包造成。当两 个红细胞或红细胞和鼓包之间的f 白j 隙小于或等于一个格子单位时，我们在它 们之间引入了润滑力。通过模拟两个对称的红细胞最后平稳地堵在鼓包之 前，验证了我们的程序的可靠性。计算结果显示，鼓包间隙6 存在两个临界 值b o o 和6 。当b c 。6 b c 时，红细胞在经过鼓包时，红细胞和红细胞之间 以及红细胞和鼓包之间可能有强的吸引作用，红细胞在经过鼓包之后有朝f f i l 管壁运动的趋势。而当d 6 1 0 时 叫罹与p o i s e u i l l e 流结果不一致，边界上速度不为零。 j 二是，有不少边界处理方法被提出希望在无滑边界上得到二级精度的结 果f 6 2 + 6 5 】。实际上这些模型只能处理直线边界。z i e g l e r 发现【5 8 】，当将 边界处_ f 流体格点和反弹格点的中间一半处( 这在后来被称为h a l f - w a y b 。n c eb a c k ) 时，可以得到二级精度的结果。半程反向弹回边界条件保留了 反向弹回边界条件便编程的优点。半程反向弹回边界条件的n o s l i p 点处于 l i n k 的中点。由此可以看出，当边界平行方格的边且处于边界l i n k 的中点 时可以实现n o s l i p 。但当边界与方格的边斜交或者边界不是直线时，边界 l i n k 的中点是阶梯状的。显然在边界上实现1 i 了n o s l i p 。 l a d d 在用晶格玻尔兹曼方法模拟粒子沉降时【1 2 】，建立了动边界条件。 将固体边界放在格点连线中间。这也就是所谓的半程反向弹回，或被称为涟 线上的反向弹回边界条件f b b l ：b o u n c eb a c ko nt h el i n k ) 。重要的是分布函 数不只是反弹，还应加上由于边界的运动而获得的动量增加的贡献： ；( i 6 ，f ) = ，：( i ，f ) 一6 w ，p ( i ，f ) 巨再6 ， ( 2 5 6 ) 其中z ( 王，r ) 是边界l i n k 上流体点碰后分布函数，z ( 磊，) 是同一l i n k 上固 体点到流体点的边界函数，玩是边界速度。用b b l 方法做粒子在流体中的 运动时由于粒子边界是阶梯状的，因此粒子的有效半径和真实半径是不_ 样 的。通常的做法是用s t o k e s 速度来测量粒子的有效半径。三维情况下的 s t o k e s 速度： k = 石赢m g ， 其中v s 是半径为r 质量为m 的小球在重力，馏作用f 在粘滞系数为u 的无界流体中下落时的稳定速度。显然这样做是比较牵强，所以后来许多人 用外推的办法建立了曲线边界条件 4 5 ，6 1 ，6 6 。 2 3 2 压力边界条件和速度边界条件 除了无滑边界条件，有时候我们要研究流体在压力差的驱动下的流动和 边界有滑移的时候流体的流动。z o u 和h e 提出【6 7 ，住反向弹回边界条 件中，在满足质量守恒的前提下，修改各个边界格点上的分撕j 函数，使之满 足一定的流体压力或速度。如图2 8 所示，对于2 t i 入口，在流动时步后分 布函数 、店和而是未知的；同理对于右出i ，在流动时步后分和函数五、 庀和一是未知的。我们的目标就是要求出在一定的压力和流速卜它们的值。 由质量和流速的定义式( 2 6 ) 、( 2 7 ) 可以得到 b 斗 l 七 l 七 ，斗 i + s + ，七f 1 七 = p 1 压力边界条件 一 、七 t 一 j 斗 s 一 1 = p u 。， 1 一 4 七 s f 1 七 b 一 = p u 。 r 2 ，5 7 ) r 2 5 8 ) r 2 5 9 ) 1 ) 在左入口，除去上下角，流动后，矗、正、五、五、兀和月已知，i 、 和 未知。设左入口的密度是p = p 。，流速度未知，竹= 0 。由( 2 5 7 ) 到( 2 5 9 1 式可以把地求出： “，：1 一盟量丛些丛丛盟。 ( 2 1 6 0 ) p l ” 6 7 48 0 u t l e t i i l l l l l l l l l l l l l l l l l l l l l l i l l l l l l l l l l l j l l l f l l 刚2 8 压力边界条件和速度边界条件示意图。 只有二二个方程要解出4 个未知量是不够的。由于 垂直于边界，因此 假设其非平衡部分满足b o u n c eb a c k 边界条件： f 、一f ；“n = f 2 一f ”j 。 代入边界处的局域平衡分粕函数，得到： 工= j + ( _ ；”一矗”) = 厶+ 6 p 。w 1 “。( 2 6 1 ) 把( 2 6 1 ) 代入( 2 5 8 ) 和( 2 5 9 ) ，分别得到： 六= 一；( 厶一) + i 1 - 3 w ：) p 。“，。 ( 2 6 2 ) = 六+ ；( 一) + ( 互1 - 3 w , ) p 。虬。 ( 2 6 3 ) 2 ) 在左入口下角，流动后，局、正、矗和力己知，_ 、五、 、石和矗 末知。由于上下是无滑边界，因此有： 。= u 。= 0 。 同理认为垂直于边界的分布函数非平衡部分满足b o u n c eb a c k 边界条 件： = 厶十( 一刚一矗圳) ， = ，4 + ( 矗一一”) 。 得到： i = f 3 + 6 p 。w l u 。， 2 2 j 4 + 6 p t n w i hv 。 因此有石砺，正气 。代入公式( 2 5 7 - 2 5 9 ) 所有未知分布函数都可以 求出： 产 ，t = = 晏砧。一( o + 一+ + + 厶+ 六十 ) 】。 左j ? 角、右上角、右下角可以用类似方法求出，只不过右边要用右边出 i a 的密度p 。代替左边入口密度p 3 ) 在右出口，除去上下角，流动后，矗、i 、正、再、石和 已知， 、 五和办未知。设右出口的密度是p = p 。，流速度毗未知， 蜥2 0 。同理 由f 25 7 ) 到( 2 5 9 ) 式可以把如求出： 圹小地逝掣。 囚此未知j 个分布函数山下列三式分别给出： = _ 一6 p 。w 】“， 2 速度边界条件 五= 六一吾( 一石) 一咕一3 w f ) 虬， 厶= 兀一圭( 一厶) 一( 1 _ 3 w t ) 成。“。 已知边界速度要求流动之后的分布函数，方法和压力边界类似。先从 ( 2 5 7 2 5 9 ) 式求得边界处和密度，再用非局域反弹求出垂直于边界的分布， 结合( 2 5 8 ) 式和( 2 5 9 ) 式则所有未知分布都可以求出来。 比如在左入口，除去上下角，流动后，矗、正、五、再、 和石已知， 小五和 未知。设矗入口的流速是度( “。，u y ) 。由( 2 5 7 ) 到( 2 5 9 ) 式可 以把密度p 求出： p ：j 一：五厶五! ! 五：五五! 。 1 一“， 用非局域反弹求出垂直于边界的分布： z = 厶+ 6 p w i “，a 结合( 2 5 8 ) 式和( 25 9 ) 式则所有未知分布都可以求出： 六= 一吉( 正一) + j 1 删，+ ( 圭一3 w ，) ， = 兀+ ；( 一 ) + j 1 彤，+ ( 1 _ 3 w i ) ，。 其它的边界情况在此不介绍，可以参考文献【6 7 】。 2 3 3曲线边界条件 以上介绍的边界条件在处理平行于格子的不动的直线边界时，是完全胜 任的。但是任我们的实际工作中往往会遇到斜交边界甚至是曲线边界。比如 研究粒子在流体中的运动、血液流中可以形变的血细胞、弹性血管等问题。 这方面有很多工作【4 5 ，6 1 ，6 5 6 9 ，这里我们只介绍f i l i p p o v a 和h i n e l 最 近提出的曲线边界条件 6 8 1 ，和m e i 等人做的一些改进 6 9 l ii 龟薯? 一 审一夺一 图2 9 曲线边界示意图。 如图2 9 所示，一条曲线边界存在于格子单位长度为斑的格点之阳j 。 一对格了，如果它们之间有一个以上的1 i n k 与曲线相交，那么固体和流体 一边的格点分别用瓦和i ，来表示。毛处的格点称为边界格点ai r 格点处碰 撞后分竹函数z ( i ，) 是已知的。对流体一边位置为i ，处的格点，在流动步 骤时我们需要知道从崮体点i 。处流到流体点i ，的分布函数z ( 毛，) 。所谓的 边界条件就是要把歹：( 毛，f ) 求出来。反向弹回边界条件只是简单地让z ( 毛，f ) 的值取成z ( i ，r ) 的值。图中 i ，= i 一贾，e 4 ，= 一i 。a i 。处的实心点是物理边界与毛和d 的连线的交点，此处流体的速度为 i 。，定义 a =0 as 1 。 ( 2 6 4 ) f i l i p p o v a 和h a n d 用一个虚构的分布函数，和i ，点的碰撞后分布函 数z ( i ，) 进行线性内插，然后再加上由于边界的运动而引起的动量增加项 一一 州 + 一一 。薛 一 一h 一。h 二一 卜巴怿一 来计算_ ，；( x b ，r ) 。 z ( i 6 ，f ) = ( 1 一z ) z ( 巧，f ) + 石1 ( 骨。，) 一6 w , p 爵舀。，( 2 6 5 ) 其中瓦是物理边界处的速度，z 是z ( 王，) 和、f ( 毛，) 之间的线性内插比 例。虚构的分布函数，”( 五，f ) 是一个类似于平衡分布的函数： 砒讲， u b s 十一耕( 2 s s ) 上式中町是i ，点流体流速，u 4 h f 是虚构的瓦点流体流速，可以外推求出： 驴等妒- d 1 弧 ( 2 6 7 ) 2 百“，+ “w 。 7 j 由三机：掣，得到： o x ，o e 。 “蛳2 + 蛾毒吣 ( 2 6 8 ) 一 i u ( 2 8 ) 式略去高阶项可以得到： 毒产惦力。刖毒坼， ( ：石，) 假设流体符合“慢流”条件： 叫( c 7 ) 1 ， 这单的上是流体的特征长度，7 1 是流体的特征时i 司，c 是格子速度。流体的 特a i 时间尺度7 1 远大于在格子尺度上的转移时间。由( 2 1 4 ) 、( 2 1 5 ) 和( 2 2 1 ) ，并略去时变音珏分，只保留空i n 变化，得到碰撞后分布函数： 才( 弓，) ：( - 7 s ，) 地枷岛昙( , t s , t ) 。 ( 2 7 0 ) 山( 28 ) 和( 2 6 5 ) 可以得到： ，1 ( 毛，f ) = 。( i ，f ) + 3 p w , i i 。( 一i ，) 。 ( 2 7 1 ) 把( 26 8 ) 、( 2 6 9 ) 、( 2 7 0 ) 和( 2 7 1 ) 式代入( 2 6 5 ) 式： z ( i 。，r ) = f 。砷( i ，) + 3 z 一3 0 一z ) ( 一1 ) 一6 a 6 p w , e , j e a ? 一蜥 。呵( 2 7 2 ) 一6 p w , e , j u 妒 在( 2 2 1 ) 式中忽略时变部分，只保留空间变化并代入( 2 1 5 ) 式，结合 ( 2 6 9 ) 式刈以得到： ，( 孑，) = ，( 膏，) 一3 咖昙a ( 2 7 3 ) 当a 趋近于1 时，应能恢复到动边界的b o u n c eb a c k 边界条件，把 ( 2 ，7 3 ) 代入( 2 5 6 ) 式得： z ( 磊，) = 。唆刁，) 一3 占细。一6 # ：7 1 4 e i e o u w j ( 2 - 7 4 ) 比较( 2 7 2 ) 、( 2 7 4 ) 式有： 3 z 一3 0 一z ) p 1 ) - 6 a = - 3 r 。 解得： z ：型。( 2 7 5 ) 存推导过程中，如果把空间导数看为一阶小量，那么上述推导达到二阶 精度。在外推虚构流速i 。，时，如果外推比例 i 2 ，则由于类似杠杆作用 误筹将被放大，a = 0 时误差达到最大而趋于无穷大。因此在a 1 2 时不宜 外推。f i l i p p o v a 提出： 咿是：!蕃1焉篙11 弼 一 2 一 m ”7 “，z = _ ，js = - 当a 。5 ，对于吉，有k i 0 5 范围内可能会出现不稳定。f 1 3 则是稳定的。如果 受模拟小r 的情况，f i l i p p o v a 曲线边界条件显然不能胜任。它只擅长于入r 的情况。 曩二 b r 1 ， n 图2 1 1f i l i p p o v a 等人的曲线边界条件的稳定范围。 后来m e i 等【6 9 提出了改进的办法。的外推公式改为下式 3 7 u b f2掣舀，+ 巩，z ：型， 。f 2 一l “ z 2 j i ， 当要时，：，由于： 一埘+ 晚暑即 ( 2 7 9 ) 由图2 9 矢口：“晒一“直+ “西一“肝= “斯一“凹= 0 ，即“晒一“口= 一( “口一“舒) ，代 入( 2 7 9 ) 式，得到： “嗡刮口地l k 言“k q 8 0 ) 把( 2 8 0 ) 式代替( 2 6 8 ) 式则( 2 7 2 ) 式化为： z ( 和) 卅”( 和) + f - 3 一( h ) ( f 咖一6 1 5 , e o 毒雄( 2 8 1 ) 一6 p w ，e 口u f p 对比( 2 7 4 ) 、( 2 8 1 ) 式即可求出z ： 2 一1 z 2 i 了。 盟然，在f 接近2 时z 趋于无穷大。同样用p o i s u i l l e 流来测试，如图 2 1 2 所示，在1 , 2 5 f 0 5 时，m e i 的曲线边界条件是稳定的。但当r 火 于1 ，2 5 时，该边界处理方法可能会发散。在小弛豫时间罩，m e i 的方法 的稳定区面积明显要比f i l i p p o v a 的方法宽。m e i 的方法和f i l i p p o v a 的 方法两者刚好互补，一个适合于小弛豫时问，另一个适合于大弛豫时问。 叻q 当 当 a2 0 5 o 0 5 0 0 0 00 10 ，20 3 0 40 5 图21 2 用故进的曲线边界条什( ；日j ) p o i s u i l l e 流进行模拟 r 的稳定范围和f i l i p p o v a 的互补。 当然，当 x ， a r c “。 图2 1 5 圆弧边界。 0 。00 51 。o 1 52 0 p ( 兀) 图2 1 6 圆弧边界上水平方向流体作用力。插图是和解析解的误筹。 4 3 2 1 0 以 童 o 4 p ( 兀) 图2 1 7 圆弧边界上垂直方向流体作用力。插图是和解析解的误差。 图2 1 5 是稳定、且静止密度为1 的流体中的一条圆弧边界。我们来 比较一下两种方法所计算得流体压力的差异。从图2 1 6 和2 1 7 中我们可 以看出，无论力的大小和分量，动量交换法都有比较明显的噪音。但是当圆 弧封闭时，误差又小了，这是对称性带来的结果。因此动量交换法模拟一般 的曲线边界也是有问题的，但是当曲线高度对称时误差又会减小，比如圆球 或者椭圆。 用动量交换法来计算流体的作用，从局部上来说是不正确的，只有大范 围平均结果是正确的。因此当所研究的问题涉及到局部的形变的时候，比如 生物膜，动量交换法原则上不能再使用。这时应该使用压力张量积分法来计 算流体对边界的作用力，虽然压力张量积分法比起动量交换法要复杂得多。 2 5润滑理论 如图2 1 8 中的方匡所示。在模拟运动边界的时候，当两个边界靠得太 近时，由于整个空间是格子化的，这时在两个边界之间可能会没有流体格子。 如此一来，将丢失流体的作用力。 图2 1 8 两个边界靠得太近时流体格子的缺失。 按照润滑理论，当两个边界靠近的时候，边界之间的流体被挤压。由于 流体不可压缩，在两个界面之间将会产生很大的压力。这个压力阻止两边界 进一步靠近。但是由于流体格子的缺失，不但压力增加不起来，反而是要减 小。进步造成界面靠近，产生恶性循环。就目前所知，解决这个问题有两 个办法。第一种是采用可变格子技术。这种办法格子不是均匀的，相当于有 p l i s a 中的自适应格子。当两个边界靠近时不断地插入新格子，以保持两个边 界之问不出现格子缺失。这种方法效率比较高，但程序复杂。如果不进行技 术处理有造成计算机内存耗尽的风险。第二种方法相对简单得多。当两界面 靠近到一定程度( 比如间隙小1 ：1 个格子) 此时可能会有流体格子缺失， 在界面引入润滑力补偿由于流体格子的缺失而损失的流体压力。l a d d 【1 2 】 是在1 9 9 4 年首次在用l b m 数值模拟粒子沉降时把润滑力引入到晶格玻 尔兹曼方法中的。 4 7 2 5 1二维圆盘之间的润滑力 y u a n 等人求出了两个等半径的二维圆盘的润滑力【7 8 】，我们将他们的 方法推广到半径不等的情况。当流体处于外力场，中时，n a v i e r - s t o k e s 方 丝0 t ，罄0 x 一古嚣0 x + ，仙等0 x ， c z s 。， 。 p 。 ： 、。 其中f 是单位质量流体所受的体积力。取质量为p 的一团流体来观察， ( 2 8 9 ) 式两边乘上p 之后，显然就是一个动量定理的方程。左边是动量的增 加，右边第一项是压力的作用，第二项是体积力的作用，第三项是粘滞力的 作用。因此有【7 9 】： 鬻。c 志v v 。蔬v u v i 。c 兰t ) t 一 粘滞力 2 i ，2 。， + 鬻骠：等u 2 0 c 。挚1 c v u ， 一= 一= 2 一一i ， 粘滞力d v 2 开 口，2 才， u 其巾，u 、l 和t + 分别是特征速度、长度和时间。加“”是无量纲纯数。 表示“j 卜比于”。r e 是雷诺数( r e y n o l d s ) 。n s 是斯托克斯数( s t o k e s ) 。 研；丝：土， s t 是斯特鲁哈利( s t r o u h a l ) 数。当雷诺数很小时，迁移惯性力可以忽略， f 2 8 9 ) 式化为： 堕o t = 一土p 竺0 x 小u 鲁o x 。 ( 2 9 0 ) ，。； 如果流动的s t 数不太大的话，雷诺数很小意味着n s 数也很小，局域惯 性力可以忽略，( 2 9 0 ) 式化为： 土p 芸o x = ，+ u 等。 ( 2 9 1 ) ，积j ( 2 9 1 ) 式就是s t o k e s 方程，满足( 2 9 1 ) 式的不可压缩流体的流动称为 s t o k e s 流动。润滑理论是以s t o k e s 流动为基础。 二受 一： j d 一k 心一_ 7 一 ) 图2 1 9s t o k e s 流中的两个二维圆盘。 如图2 1 9 所示，s t o k e s 流中两个靠得很近的二维圆盘。几何参数和 速度如图中所示。队为压力沿y 方向无变化，s t o k e s 方程化为： _ d p ：娶，( 2 睨) 出 a y 。 其- ，速度分解为( “，v ) ，“是水平分量，v 是垂直分量。边界条件为： 归，、挈。 ( 2 9 3 ) v2 v a 刍y 2 一“h吼心r + + ，、【 = 定义r a y l e i 曲耗散函数 其中能量耗散率 其中形变率张量为 相应的耗散力可以求出 孵；三胁， 1j d = 2 , u d ：d ， d ；v 开。 f ：一塑。 d ” f 2 9 4 ) r 2 9 5 ) r 2 9 6 ) 由于v “，忽略掉速度的垂直分量，而且娑 0 ，5 ： d o c ) 如果没有一个方向满足条件，则o 点的分布函数就直接取最近流体 格点的分布函数。如果有多个方向满足条件则取平均。 实践证明，这样可以减小误差。事实上我们发现当雷诺数很小的时候， 直接取最近流体格点的分布函数作为边界上的分布函数也可以。模拟中圆柱 的半径我们取1 3 个格子单位，弛豫时间f = 0 6 。 0 。6 0 0 5 5 0 s o 0 4 5 j0 4 0 0 3 5 0 3 0 0 2 5 o 。2 0 - - - 。 - 七1 = _ ! 避曳j 02468 1 0 y ，l 图4 2 不同雷诺数，圆柱质心轨迹。实线是有限元计算的结果。 图4 2 是不同雷诺数下，圆柱的质心轨迹。可见我们的计算和有限元的 结果符合得很好。雷诺数大的时候出现过冲现象。 迎 3 ，h 望 e u 一 t ( s ) t ( s ) t ( s j t ( s ) 例4 3不同雷诺数f ，圆柱质心水平方向速度随时间的变化。 图4 3 是不同雷诺数下，圆柱质心水平方向速度随时问的变化。用了 压力张量积分法，水乎方向速度与有限元法精确吻合，解决了图3 6 中动 量交换法计算得的水平方向速度与有限元法符合得不好的问题。从图2 4 、 2 1 6 和2 1 7 我们可以看出，由于速度的离散化，会造成一些方向动量的丢 失或者多出一些方向的动量来。如果边界是直线，长直线的平均结果是正确 的。如果边界是曲线，在格子上又不对称的话，多出的和丢失的不能互补， 则动量的计算是有较大的误差。张量积法是把离散了的速度再连续化，所以 流体作用在边界上的作用力的计算精度将会提高。 芝 e u t ( s ) t ( s ) 鼍 誊 t ( s ) t ( s ) 图44 不同雷诺数下，圆卡_ ：质心垂直方向速度随时间的变化。 图4 _ 4 是不同雷诺数下，圆柱质心垂直方向速度随时问的变化，也比 动量交换法的结果好。 一，o 一 一s ，e u ) 芝 可 爱 8 t l s ) t ( s ) 图4 5 不同雷诺数下 t ( s ) t ( s ) 圆柱角速度随时间的变化。 图4 5 是不同雷诺数下，圆柱角速度随时问的变化。角速度也比动量 交换法的符合得好。总之压力张量积分法改进了计算的精度。 4 2单个悬浮