微通道中液滴的耗散粒子动力学模拟

　第39卷第11期　2005年11月

上海交通大学学报

JOU

RNAL O F SHAN GHA I J I AO TON G UN I V ER S IT Y

V o l . 39N o. 11　

N ov . 2005　

　　文章编号:100622467(2005) 1121833205

微通道中液滴的耗散粒子动力学模拟

陈　硕1, 　赵　钧1, 　王　丹1, 　刘明安2, 　范西俊3

(1. 上海交通大学航空航天工程系, 上海200030; 2. 空军工程大学电讯工程学院, 西安710077;

3. Cooperative R esearch Cen tre fo r Po lym ers , 32B u siness Park D rive ,

N o tting H ill , V ic . 3168, A u stralia )

摘　要:利用耗散粒子动力学(D PD ) 对三维微通道中的液滴动力学进行了研究. 通过提高两种流体之间的保守力系数获得不互溶的两种流体, 并分析了相应的F lo ry 2H uggin s 模型中的ς参数.

. 计算采用消息传递并行计算技术D PD 计算所获得的表面张力与Groo t 等给出的理论公式相符合

(M P I ) , 计算过程模拟了高分子液滴穿过突扩突缩通道时的形态变化. , 不同尺寸的高分子液滴穿过微通道所需要的时间不同, 升, 但上升的趋势逐渐减缓.

关键词:耗散粒子动力学; 液滴; 微通道; 中图分类号:O 35　　　文献标识码:A

The D P D ynam ics S i m ula tion of

in M ic rocha nne l

S huo , 　ZH A O J un , 　W A N D an , 　L IU M ing 2an , 　FA N X i 2jun

1

1

1

2

3

(1. D ep t . of A eronau tic &A stronau tic Eng . , Shanghai J iao tong U n iv . , Shanghai 200030, Ch ina ; 2. T elecomm un icati on Eng . In st . , A ir Fo rce Eng . U n iv . , X i’an710077; 3. Cooperative R esearch Cen tre fo r Po lym ers , 32B u siness Park D rive , N o tting H ill , V ic . 3168, A u stralia )

A bs tra c t :T he flow dynam ics of drop s in a 3D m icrochannel w ere investigated u sing the dissi pative p article

dynam ics (D PD ) m ethod . T he N ew ton ian drop is m ade up of si m p le D PD particles , and the po lym er drop is m ade up of m any fin ite ex ten sib le non 2linear elastic (FEN E ) bead sp ring chain s , w ith each chain con sists of 16beads . T he so lid su rface of the m icrochannel is m ade up of “frozen ”si m p le D PD particles . T he b inary i m m iscib le flu ids w ere si m u lated by increasing the rep u lsi on coefficien t betw een differen t flu id p articles , the ςp aram eters in F lo ry 2H uggin s m odels w ere also analyzed acco rdingly .

T he value of

in terfacial ten si on com pu ted by D PD is in agreem en t w ith Groo t and W arren’stheo retical p redicti on . T he confo r m ati on evo lu ti on of po lym er drop p assing th rough con tracti on 2diffu si on m icrochannel w as also described . T he p arallel com pu tati on M P I techno logy w as i m p lem en ted in the com pu tati on p rocess . T he com p u tati onal resu lts show that the ti m e needed to pass th rough the slit m icrochannel is differen t fo r differen t drop sizes , and the ti m e in itially increases rap idly m ono ton ically w ith increasing of the rati o of M H (diam eter of drop heigh t of m icrochannel ) , and then gradually increases slow ly .

Ke y w o rds :dissi p ative particle dynam ics (D PD ) ; drop ; m icrochannel ; parallel com pu tati on

收稿日期:2005204225

基金项目:国家自然科学基金资助(50576051) , 上海交通大学青年教师校内科研启动基金资助项目

作者简介:陈　硕(19692) , 男, 四川泸州人, 副教授, 主要研究方向为微尺度流动与传热数值模拟. 电话(T el . ) :[1**********]4;

E 2m ail :chenshuo @sjtu . edu . cn .

1834

　　　　上　海　交　通　大　学　学　报

第39卷　

　　耗散粒子动力学(D PD ) 方法是一种介观尺度的数值模拟方法[1], 该方法适合于研究复杂流体的动力学行为. 在D PD 系统中, 基本的单元是一些离散的被称为粒子的动量载体, 这些所谓的粒子具有粗粒化的概念, 每颗粒子的运动代表的是大量分子(或原子) 的集体行为, 而非分子动力学模拟中跟踪单个分子(或原子) 的运动行为. 相对于分子动力学中分子(或原子) 之间的势能来说, 耗散粒子动力学中粒子之间的相互作用要“软”得多, 因此, 耗散粒子动力学所模拟的时间和空间尺度要比分子动力学大得多. D PD 还可以理解为宏观的微分流动控制方程在小尺度上的随机描述[2], 从这个意义上说, D PD 方法是连接微观分子动力学方法和宏观流体力学方法的一座桥梁, 是一种真正的介观尺度的模拟技术.

液滴动力学是多相流体系统研究的一个重要方面. 液滴的动力学特性影响着大量的天然与人工过程, 许多与工业、生物有关的系统, 如乳状液、聚合物共混物以及生物流体等都可以被认为是液滴的悬浮物. 形、破碎以及聚合等行为[3～6]. . , 研究结果将为与液滴分离有关的M E M S 器件的设计提供一定的依据.

式中的求和符号对一定的截断半径r c 内所有的不

包括粒子i 的其他粒子进行求和(如果粒子之间的距离r ij >r c , 则粒子之间的相互作用为零) .

保守力F C, ij 为作用在2个粒子质心连线方向上的粒子之间的排斥力, 可表示为

F C, ij =

a ij 1-δ

r ij r c

r ij

r ij

(3)

式中:a ij 为保守力系数, 表示粒子i 与粒子j 之间最δij =大的排斥力值; r ij =r i -r j , r ij = r ij , r

粒子j 作用在粒子i 上的耗散力(或者阻力)

F D , ij 可以表达为

δδ(4) F D , ij =-Χw D (r ij ) (r ij v ij ) r ij 式中:w D 为与粒子间距离r ij , 当r ij >

r c 时, w D =0; v ij , v ij =v i -v j ; Χ, Χ, . 系统动能的降低将由随机力F R , ij 引起的粒子的随机运动进行补偿, 随机力可表示为

δF R , ij =Ρw R (r ij ) Νij r ij

(5)

式中:w R 为与粒子之间距离r ij 有关的权函数, 称为

随机力权函数, 当r ij >r c 时w R =0; Νij 为随机变量, 其平均值为零且方差为∃t -1, ∃t 为时间步长; Ρ为随机力系数. 与保守力一样, 随机力和耗散力也作用在粒子质心的连线方向.

[7]ζ(5) 中的2个E sp an o l 等的研究表明, 式(4) 、

权函数相互关联, 两者只能任取一个, 其相互关系可表达为

2

(6) w D (r ij ) =w R (r ij ) 耗散力系数Χ和随机力系数Ρ也相互关联, 两者也只能任取一个:

Ρ2=2Χk B T

(7)

1　耗散粒子动力学

1. 1　数值方法

基于牛顿运动方程, D PD 系统中粒子的位置和

速度随时间的变化规律为

(1) =v i , 　=f i +F e

d t d t

式中:r i 和v i 为粒子的位置和速度矢量; f i 为所有其他粒子(不包括粒子i 本身) 作用在粒子i 上的力; F e 为用于模拟高分子链的弹性力或为所受到的外力等; t 为时间. D PD 系统中单位质量为单个粒子的质量, 因此作用在每个粒子上的力的大小等于粒子的加速度值. 粒子之间的动力相互作用由耗散力和随机力组成, 这两部分作用相互补偿使得系统的平均动能为一恒定值.

耗散力F D , ij 和随机力F R , ij f i 包括保守力F C, ij 、三部分, 其中每一部分都表示粒子之间的两两相互

作用:

f i =

式中, k B T 为系统的Bo ltz m ann 温度, 与系统的涨落-耗散理论相类似[8]. 将k B T 作为能量的单位, 则有:

(8) Ρ2=2Χ　　采用的权函数表达式为

w D (r ij ) =w (r ij ) =

2

R

1-0

r c

r ij

r ij ≥r c

(9)

∑

j ≠i

(F C, ij +F D , ij +F R , ij ) (2)

这种形式的权函数能够在粒子之间产生较强的耗散力. 在本文中, 截断半径r c 取为单位长度, 即r c =1.

　第11期

陈　硕, 等:微通道中液滴的耗散粒子动力学模拟

1835

1. 2　两相不互溶性流体的模拟

在D PD 中, 两种流体之间的互溶性主要由两种流体之间的保守力系数a ij 决定. 为了模拟两种不互溶的流体, 按照Ro thm an 2Keller 的方法[9], 引入一种新的被称为“颜色”的变量,

a ij =

a 0　i 和j 为同一种颜色粒子a 1　i 和j 为不同颜色粒子

应力张量的计算可以通过Irving 2K irkw ood 方

法[11]得到:　S ΑΒ=-V

(10)

当不同颜色的两种粒子相互作用时, 如果改变它们之间的相互作用的保守力系数, 则改变了它们之间的排斥力. 如果a 1≈a 0, 则两种流体完全互溶, 如果

. a 1µa 0, 则两种流体不互溶

[10]

Groo t 等给出了保守力系数a ij 和F lo ry 2

式中:m i 为粒子的质量(本文中m i =1) ; N p 为粒子的数量; u i Α和u i Β为粒子i 的本动速度在Α、Β方向的

λΑ(x ) , v λΑ(x ) 为在位置x 处流场分量, 例如u i Α=v i Α-v

速度Α方向的分量; F ij Β为由粒子j 施加在粒子i 上

的力Β方向的分量; V 为体积; r ij Α为粒子i 和j 之间距离矢量在Α方向的分量; 〈…〉为求平均. 式(12) 右边的第1项为D PD 粒子的动量传递对应力的贡献, 第2项为粒子之间作用力的贡献. 压力可表达为

(13) p =-tr S

3

式中:S 为应力张量, 由式(12. 1. 5N p

i

N p N p

ij Α

∑m

i

u i Αu i Β+

∑∑r

i

j >1

F ij Β

〉

(12)

. 如果A 和B 为两种H uggin s 模型中ς参数的关系

不互溶性流体时, 则ς值为正, 否则ς值为负. 当ς

值为正且超出一定的临界值时(即ς>ςcrit ) , 则会出现两相分离, 粘性流体中悬浮的液滴就属于这种情况. ςcrit 可由Groo t 等[10]给定的公式进行计算. 在本文中相应的参数选择使得两相分离发生, 也即ς>

[10]

ςcrit . Groo t 等给出了简单流体ς种密度(Θ=3, 5) 下ς∃a ij 的函数, 其中∃a ij =10. , ∃a ij 0. ) a ij 　　Θ=3ς=

(0. 0. 002) ∃a ij 　　Θ=5

流体B #AB 可以通过

[]

得到:

AB =

p z z -

(p x x +p y y ) d z 2

(14)

(11)

在本文的计算中, D PD 系统的密度值为4, 因此可以利用式(11) 估计Θ=

4条件下ς值的范围. 1. 3　高分子液滴的构成

D PD 方法中, 由牛顿流体所组成的液滴可以由

式中, p x x 、. 本文p y y 、p z z 为压力张量的3个对角分量

假定界面与xOy 平面平行, 如图2所示. 表面张力可以进一步表达为

(p x x 〉) (15) #AB =L 〈p z z 〉-〈+

〈p y y 〉2

式中, L 为所模拟区域的高度. Groo t 等[10]给出了计算简单流体表面张力的理论公式:

0. 26±0. 01

#AB =(0. 75±0. 02) Θk B T r c ς[

1-2

(2. 36±0. 02) ς]3

简单的D PD 粒子构成, 只需提高液滴与周围流体之

间的保守力系数即可. 本文采用多根有限拉伸非线性弹性珠簧链(FEN E 链) 相互纠缠来构造出高分子液滴[5], 如图1所示. 其中, 单根的FEN E 链包含16个珠子, 液滴中10%的粒子为简单粒子, 而其余90%的粒子由多根FEN E 链组成. 对于松弛状态下半径为5的液滴, 如果液滴密度为4, 则液滴内部包含2096个粒子.

(16)

图2　表面张力计算示意图

F ig . 2　Sketch diagram of interfacial tensi on computati on

2　针对D PD 的并行计算技术

耗散粒子动力学本质上属于离散粒子动力学方法, 它要求巨大的计算量才能展现真实流体的物理特性, 一个有效的解决办法是采用高性能并行计算技术. 针对耗散粒子动力学所描述的系统, 可以采用消息传递M P I 并行计算技术来实现对较大系统的

图1　由多根FEN E 链所构成的高分子液滴F ig . 1　A po lym er drop m ade up of FEN E chains

1. 4　应力计算

1836

　　　　上　海　交　通　大　学　学　报

第39卷　

真实模拟. 和被分配到各个微处理器上的粒子数量相比, 参与通信的粒子数量是很小的一部分, 因此可以保持较高的并行效率.

3　计算结果分析

各参量均由D PD 单位进行无量纲化处理[5]. 在本研究中, 采用了14个微处理器实施并行计算, 其效率维持在90%以上. 当∃a ij =37. 5时, ς值为10. 65～25. 91. 本文中D PD 的计算方法采用Groo t

[10]

等给出的velocity 2V erlet 方法. 为了计算表面张力, 将计算区域设置成为40×10×30, 流体A 为组成液滴的粒子, 流体B 为组成周围流体的粒子(见图2) . 牛顿流体液滴和高分子液滴表面张力的计算结果如表1所示.

表1　表面张力D PD 计算值

Tab . 1　Co m put ational value of i n terfac i al ten sion by D PD

a 1 a 0

将开始运动并被挤入比其本身直径要小的微通道内

且穿过该狭窄通道, D PD 的模拟结果如图4所示. 图4中只画出了组成液滴的粒子, 组成周围流体的粒子未画出.

高分子液滴

0. 0910. 8442. 5853. 5505. . 牛顿流体液滴

0. 0160. 0211. 540760. 000

1. 01. 21. 51. 83. 05. 0

图4　高分子液滴穿过狭缝通道的D PD 模拟结果

F ig . 4　D PD si m ulati on

results of po lym er drop

passing th rough slit channel

, ∃a ij =37. 5时的　　根据式(:

Θ=3时#AB =3. 81; Θ=5时#AB =6. 06. 由此插值得到Θ=4时

#AB =4. 94, 与本文D PD 模拟中∃a ij =37. 5、a 1 a 0=

3. 0、=4时, #AB =5. 25相对应. 由表1可见, 表面Θ

　　图5所示为不同尺寸的高分子液滴在穿过微通

道时所需要的时间t 和M . t 为穿过H 之间的关系Q =20所需要的时间, H =6. 随着M H 的增加, t 单调上升, 但当M . H >1. 3时, t 变化不大

张力随着a 1 a 0值的增大而增大; 在同等条件下, 高

分子液滴的表面张力大于牛顿流体液滴的表面张力.

由此进一步计算了高分子液滴通过微小通道时的情况, 如图3所示, 并对x 、y 方向实施了周期性边界条件. 图中:M 为松弛状态下液滴的直径; H 为微通道的高度; Q 为微通道测试段的长度.

图5　M H 与t 的关系

F ig . 5　T he relati on betw een M H and t

4　结　论

图3　液滴穿过微通道示意图

F ig . 3　Sketch diagram

of po lym er drop passing

th rough m icrochannel

(1) 耗散粒子动力学可以十分方便地模拟液滴

悬浮物在微通道中的运动, 表明耗散粒子动力学在处理界面流动方面有明显的优势.

(2) 利用耗散粒子动力学计算得到的液滴表面张力与Groo t 等所给出的理论计算公式相符合; 在

　　微通道中固体壁面采用冻结的简单D PD 粒子构成, 对如图3所示的系统施加一定的外力, 则液滴

　第11期

陈　硕, 等:

微通道中液滴的耗散粒子动力学模拟

1837

同等条件下, 高分子液滴的表面张力比牛顿流体液滴的表面张力大.

(3) 不同尺寸的高分子液滴在穿过微通道时所需要的时间t 不同; 随着M H 的增加, t 值单调上升, 但上升的趋势逐渐减缓. 参考文献:

[1]　Hoogerbrugge P J , Koel m an J M V A . Si m ulating

m icro scop ic

hydrodynam ic

phenom ena

w ith

dissi pative particle

dynam ics [J ].

Europhysics

[6]　L iu Y , Yang X Z , Yang M J , et al . M eso scale

si m ulati on on the shape evo luti on of po lym er drop and initial geom etry influence [J ]. Poly m er , 2004,

45(4) :6985-6991.

ζo l P , W arren P B . [7]　E span

dissi pative particle

Statistical m echanics of

Europhysics

dynam ics [J ].

L etters , 1995, 30(4) :191-196. [8]　H uilgo l R R , Phan 2T h ien N .

viscoelasticity :

General

F luid m echanics of

constitutive

p rinci p les ,

modelling , analytical and num erical techniques [M ]. Am sterdam :E lsevier , 1997.

[9]　N ovik K E , Coveney P V . Sp inodal decompo siti on of

off 2critical quenches w ith a viscous phase using dissi pative particle dynam ics in tw o and th ree spatial di m ensi ons [J ]. Physical Rev iew E , 2000, 61(1) :435-448.

[10]　Groo t R D , W .

m ic si m [J ].

D issi pative particle Journal of Che m ical

ics :een atom istic and 11) :4423-4435.

], Phan 2T h ien N , Yong N T , et al . M o lecuar

dynam ics si m ulati on of a liquid in a comp lex nano channel flow [J ]. Physics of Fluids , 2002, 14(3) :1146-1153.

[12]　Irving J H , K irkw ood J G . T he statisticalm echanical

. I theo ry of transpo rt p rocesses V . T he equati ons of hydrodynam ics [J ].

Journal of Che m ical Physics ,

1950, 18(6) :817-829.

L etters , 1992, 19(3) :155-160.

ζo l P , Serrano M , Zuniga I [2]　E span . Coarse graining of

a fluids and its relati on w ith dissi pative particle dynam ics and s moo thed particle dynam ics [J ]. I n ternational Journal of M odern Physics C , 1997, 8(4) :899-908.

[3]　Jones J L , L alM , R uddock J N , et al . D ynam ics of

a drop at a liquid so lid interface in si m p le shear fields :A m eso scop ic si m ulati on study [J ]. Faraday D iscussion s , 1999, 112:129-142.

[4]　C lark A T , L alM , R uddock J N , et al . M eso si m ulati on of drop s in [J ]. Lang m uir , 15:.

[5]　Chen S , X J , al . D issi pative

particle si m on of po lym er drop s in a peri odic [J ]. Journal of Non -Newton i an Fluid M echan ics , 2004, 118(1) :65-81.

下期发表论文摘要预报

基于有源辅助的被动跟踪系统研究

李安平, 　敬忠良, 　胡士强

(上海交通大学航空航天信息与控制研究所, 上海200030)

摘　要:提出了一种有源雷达间歇辅助红外工作的雷达 红外协同跟踪系统, 该系统通过减少雷达开机时间降低被敌方ES M 锁定的可能性. 为了改善纯角度跟踪的滤波发散问题, 一方面, 在滤波处理时采用了非线性逼近能力更强的U nscented 卡尔曼滤波算法; 另一方面, 充分利用雷达、红外同时开机时的量测信息, 构造出一组时间多项式, 在雷达关机期间, 利用该组时间多项式估计目标的运动状态, 辅助红外传感器进行跟踪. 计算机仿真结果表明, 该方法不仅能保证跟踪精度, 同时也能有效地降低雷达被敌方ES M 锁定的可能性.