适用复杂几何壁面的耗散粒子动力学边界条件*

林晨森 陈硕† 肖兰兰

1)(同济大学航空航天与力学学院,上海 200092)

2)(上海工程技术大学机械与汽车工程学院,上海 201620)

耗散粒子动力学(DPD)是一种针对介观流体的高效的粒子模拟方法,经过二十多年发展已经在诸如聚合物、红细胞、液滴浸润性等方面有了很多研究应用.但是因为其边界处理手段的不完善,耗散粒子动力学模拟仍局限于相对简单的几何边界问题中.本文提出一种能自适应各种复杂几何边界的处理方法,并能同时满足三大边界要求：流体粒子不穿透壁面、边界处速度无滑移、边界处密度和温度波动小.具体地,通过给每个壁面粒子赋予一个新的矢量属性—局部壁面法向量,该属性通过加权计算周围壁面粒子的位置得到；然后通过定义周围固体占比概念,仅提取固体壁面的表层粒子参与模拟计算,减少了模拟中无效的粒子；最后在运行中,实时计算每个流体粒子周围固体粒子占比,判断是否进入固体壁面内,如果进入则修正速度和位置.我们将这种方法应用于Poiseuille流动,验证了该方法符合各项要求,随后还在复杂血管网络和结构化固体壁面上展示了该边界处理方法的应用.这种方法使得DPD模拟不再局限于简单函数描述的壁面曲线,而是可以直接从各种设计图纸和实验扫描影像中提取壁面,极大地拓展了DPD的应用范围.

1 引 言

耗散粒子动力学(DPD)由Hoogerbrugge和Koelman[1]在 1992年提出构想,随后由 Español和Warren[2]做出重要的系统论证并搭建了理论框架,经过二十几年的发展,至今已经成为介观尺度模拟流体动力学的最主要的方法之一,被称为“弥补宏观尺度和介观尺度之间空白的方法”[3].它的基本思想是一些离散的被称为粒子的动量载体在连续的空间和离散的时间上运动,这些粒子代表的是一个小区域内的大量分子的集体行为,这样就可以在粗粒化粒子尺度上进行研究,从而忽略更小尺度的结构和运动状况.经过人们的不断探索和尝试,DPD在以下应用场景中展现出独特的模拟优势: 1)聚合物,采用珠簧链模型,可以建立和Flory-Huggins理论之间的对应关系,研究微相分离体系和自组装现象,可用于药物输运、介孔材料制备等研究; 2)红细胞,DPD红细胞模型可以准确对应真实红细胞的力学特性,已经再现了单个红细胞穿过狭缝、多个红细胞聚集成串、红细胞在血管狭窄处堵塞流动等现象[4,5]; 3)液滴浸润,对DPD势能稍加修改就可以模拟气液共存的系统[6],对液滴在各种化学和几何异质表面上运动行为的研究为预防表面结冰和微流控芯片中的液滴控制带来新的启示[7,8].此外,随着近年来计算机图形处理器(GPU)作为通用计算加速器的蓬勃发展,DPD已经实现了在GPU上的加速运行,达到了比CPU快几十倍的运算速度[9],更快的计算时间意味着可以计算更大的模拟体系,这使得该方法的应用边界进一步扩大,成为愈发有力的模拟研究工具.

在近年来的DPD研究中,边界条件的设置一直是研究者们关心的热点,各种相应的处理方法被提出以应对特定的问题.发展至今,DPD模拟中对边界的处理大致分为两类.第一类是采用周期性边界条件,即不存在传统意义的固体边界,一般应用于不受限流体.在此基础上进行改进,还进化出模拟双Poiseuille流用以测量流体黏性的边界方法[10]和模拟库特流的Lees-Edwards边界方法[11]等.第二类方法是构造固体的边界,也即将问题变成受限流体问题.以上第一类周期性边界条件虽然能较好地处理某些特定的问题,但应用范围很窄,大多数流动问题是在管流内或者受到某些限制的,例如在微通道内的大分子流动[12]、多孔介质中的多相流[13]、结构化表面上的液滴行为[14]等,都需要用到有形的固体边界.固体边界条件还可细分为弹性壁面和非弹性壁面,弹性壁面指可以发生弹性变形的壁面,例如血管壁等,因为技术原因,DPD目前还不能很好地处理这种移动的弹性的固-液界面,也鲜见相关弹性壁面的文献.相比之下非弹性壁面更为简单,在模型中较易实现,是采用最多的边界处理方法.在以往的DPD研究中,非弹性壁面往往还被进一步简化成没有结构的平板或简单的几何形状,究其原因主要是缺乏对复杂几何边界的处理方法.这使得DPD方法的应用受限,在模拟进一步的复杂问题,例如红细胞在不规则血管网络中的流动现象、液滴在微结构表面的运动时无法准确地处理这样的边界.一种对复杂几何边界自适应的边界处理方法将提升DPD对微流体问题的研究能力.

2 耗散粒子动力学的边界条件

2.1 耗散粒子动力学模型

DPD模型中包含一系列的质点粒子,这些粒子在无网格的空间上运动,彼此之间根据三种力互相碰撞: 根据势能推导出来的保守力降低粒子之间切向速度的耗散力粒子连线方向上的随机力后两种力可以看作是配对的布朗阻尼器,和郎之万动力学或者布朗动力学不一样,DPD的这对力是动量守恒的.对于粗粒化的DPD系统,布朗阻尼器是在粒子之间体现出黏性力的同时还要体现热噪音的最简单模型.因为动量是守恒的,所以当尺度达到一定量级后模型中流体的行为就完全符合宏观的流体力学了.

在DPD中,粒子的运动符合牛顿第二定律,保守力、耗散力和随机力的计算公式如下:

其中 rij=|ri-rj| 是粒子i和粒子j之间的距离;eij是连接粒子i和粒子j的单位向量; vij是粒子i和粒子j之间的相对速度; θij是高斯白噪音,同时具有对称性,即θij=θji,以此保证整体的动量守恒; γ和σ分别是耗散力参数和随机力参数,它们之间满足耗散涨落定理:

其中kB是玻尔兹曼常数,T是平衡温度.简单地根据距离的衰变函数来定义保守力的权重函数,

耗散力和随机力的权重函数采用以下形式:

其中在经典DPD中s=2,s也可以取其他值用来调节流体的黏度.

2.2 边界条件的要求

根据长期的实践研究,人们总结出优秀的固体边界条件至少要满足:

1)流体粒子不能穿透壁面;

2)流体在固体壁面处无滑移;

3)固体壁面不能影响周围流体的物理性质,如固体壁面附近的温度、密度要和流体内部一致,不能波动太大.

一种简单并且被广泛应用的方法是在壁面的位置排布冻结的粒子,通过冻结粒子和流体粒子的作用来反映固体对流体的力.这种方法和MD中处理固体壁面的方法很像,但在DPD中会面临新的问题: 因为粒子之间的作用势相对较软流体粒子很容易穿透固体壁面粒子从而进入固体内部.为了解决这一明显的错误,研究者们提出了两种方法:1)提高固-液排斥作用力; 2)设置虚拟反弹边界.提高固-液排斥力一般通过提高固体粒子密度或者增大壁面粒子与流体粒子间的排斥力系数aij来实现,操作起来非常简单,对于防止穿透的效果也非常好,但人为设置的强烈壁面排斥力虽然防止了流体的穿透,也严重影响了固体壁面附近流体的物理性质,例如温度和密度,进一步影响边界附近的流动行为使模拟结果偏离正确结果.设置虚拟反弹条件则不改变固-液间的作用强度,只对模拟过程中穿透进入固体内部的流体粒子进行位置和速度的修正,主流的处理方法有以下三种: 镜像反射、原路折回反射和麦克斯韦反射.通过反复比较和实践,这三种方法均存在明显的缺点[15],且当壁面的形状较复杂时,难以界定粒子是否穿透进入固体壁面内部.

为了解决壁面附近液滴粒子的均匀性问题,Xu和Wang[16]赋予壁面粒子相对流体粒子的虚拟速度,用于计算对流体粒子的耗散力,从而增加了流体粒子的耗散阻力实现了壁面的无滑移边界条件,并得到平滑的壁面附近温度密度分布,但仍需要依靠人为指定的虚拟边界来界定哪些流体粒子已穿透需反弹,并不能自动适应复杂壁面结构.为了解决DPD对复杂几何壁面的适应性问题,Mehboudi等[17]用三角形单元描述边界形状,这样可以容易地从各种微机电系统设计图中获得边界形状,并经测试获得了很好的边界效果,但因为需要在粒子系统中引入网格处理的边界方法,所以在程序简易型和计算量上仍存在推广的壁垒.刘谋斌和常建忠[18]将整体计算域用规则背景网格覆盖,分为流体区域网格和固体障碍区域网格,固体障碍区域网格根据周围网格的种类来计算法向量,实现了对诸如多孔介质等问题的边界处理,但仍存在密度波动,并且用规则网格背景对复杂几何的描述仍需要近似.Li等[19]提出了一种模拟运行时动态检测粒子是否穿透壁面并动态计算周围壁面法向量的方法,这种方法在复杂微流控芯片内通道内展现出很好的适应性,并对旋转双筒等动态壁面的情况也能适应.对于固定和移动壁面,这种方法不加区分,运行时均需每一步都重新计算壁面法向量.

3 新型的适用复杂形状的边界条件

我们提出一种新型的可以应用在任何复杂几何形状上的边界条件LWNM(local wall normal method).具体地,新定义一个粒子的矢量属性——固体壁面局部法向量(LWN),并将其记录在每个固体粒子上,例如壁面粒子i的局部壁面法向量记为lwni.如果壁面是可移动的,在每个时间步计算固液之间的作用力前,需要先更新每个固体粒子的法向量; 如果固体壁面是静止的,则无需更新,只需在模拟初始化时计算一次即可,额外的计算开销接近于0.

图1是计算壁面粒子i的法向量lwni的示意图,其中S代表固体壁面区域,rcw是计算壁面局部法向量时的截断半径,可以不同于粒子之间作用力的截断半径rc; 固体粒子j是固体粒子i的rcw范围内的其他固体粒子,rij是粒子j到粒子i的距离,rij是粒子j到粒子i的向量.

图1 计算计算壁面粒子i的法向量lwniFig.1.Computing the local wall normal vector lwni of wall particle i.

计算固体粒子i的局部壁面法向量的方向向量lwnti:

将lwnti的长度进行标准化后得到壁面粒子i的壁面局部法向量lwni.:

为了判断流体粒子是否已经穿透壁面进入固体区域,为每个流体粒子新增一个额外属性φ,定义φ为流体粒子的固体体积分数,通过其周围的固体壁面粒子的位置计算:

其中ρw是固体壁面的平均粒子密度,是一个三维的权重函数,对rc范围内的空间进行权重的处理,这个权重函数也被应用在多体耗散粒子动力学(MDPD)的保守力计算中.如图2,当粒子完全浸入到固体壁面中时,粒子截断半径内被固体粒子充满,根据(10)式计算得φ应为1,实际中因为壁面粒子的密度波动,所得φ在1左右波动(壁面粒子密度越高,波动范围越小,密度越小,波动范围越大); 当粒子远离固体壁面时,粒子截断半径内没有固体粒子,φ为0; 当流体粒子位于理想固-液的分界面上时,半径rc球内的一半被固体占据,另一半被流体占据,根据(10)式得φ应为0.5,如果φ ＞ 0.5即对应流体粒子更多地进入固体,有更多的固体粒子在截断半径内,如果φ ＜ 0.5即对应流体粒子更远离固体,有更少的固体粒子在截断半径内.因此以φ=0.5为分界线,判断流体粒子是否穿透固体壁面.

图2 流体粒子固体体积分数(φ)的四种情况(φ=0远离壁面,0 ＜ φ ＜ 0.5 未穿透壁面,0.5 ＜ φ ＜ 1.0 已穿透壁面,φ=1.0完全浸入壁面)Fig.2.Four scenarios for solidfraction (φ): φ=0 particle away from wall; 0 ＜ φ ＜ 0.5 particle near the wall; 0.5 ＜ φ＜ 1.0 particle penetrating the wall; φ=1.0 particle submerged in wall.

采取预估-修正的策略来防止流体粒子穿透壁面.在每一个时间步,预估粒子下一时刻的位置并计算φ,当φ ＞ 0.5,即表明粒子下一个时间步会穿透壁面,需要对流体粒子速度进行修正,修正后的新速度为

其中U和a是壁面局部的速度和加速度.如果模拟问题中壁面静止,(11)式可以简化成

其中en是流体粒子对应的固体壁面局部法向量.如图3所示,在流体粒子的截断半径内,有若干带局部法向量的壁面粒子,en取其加权平均后归一化的向量,权重函数用的是MDPD保守力中采用的三维权重函数:

图3 选用一定范围内固体粒子的lwni计算修正流体粒子时的法向量enFig.3.Computing en from nearby lwni for correcting penetrating fluid particle.

在对即将穿透固体壁面的流体粒子进行速度修正时,(11)式的本质是原路返回反射方法,原路返回反射方法对壁面附近的密度和温度波动影响较小,但是并不能很好地满足无滑移边界条件.为了改进这种方法,研究者们对此提出了很多修改方法[19—23].本文采用文献[19]提出的办法对固体和流体粒子之间的耗散力系数进行修改.假设流体粒子距离固体壁面的距离为h,则修改后的耗散力系数为

图4是一个不规则流体通道LWNM边界条件实施的流程图.操作步骤如下: 1)提取固体区域的边界层粒子,如图4(b),为了简化计算,只有固体边界层粒子会被赋予局部壁面法向量,固体内部区域的粒子在LWN创建完成后可以删除,只留下壳层固体粒子,这样可以进一步减少系统总粒子数,加快计算速度; 2)通过(8)和(9)式计算壳层固体粒子的LWN,每个固体粒子都有独立的局部壁面法向量,如图4(c),红色箭头即为局部壁面法向量; 3)在模拟过程中,通过(10)式预判流体粒子下一步是否会进入固体壁面内,如果是则通过(11)式和(15)式修正粒子速度,模拟时的快照如图4(d).

图4 LWNM边界条件的实施过程图 (a)固体粒子构造壁面; (b)提取表面层固体粒子; (c)计算LWN; (d)模拟时效果Fig.4.Workflow of LWNM: (a)Constructing wall with frozen particles; (b)identifying surface wall particles; (c)computing LWN; (d)snapshot during simulation.

LWNM不仅适合CPU计算,也同样适合GPU计算,无论是LWN的计算还是之后φ的计算,都只需要截断半径内邻居粒子的信息,而这些信息是标准DPD模型中计算力都已经计算过的,无论CPU计算还是GPU计算都可以在现有的模型上经过少量的改动实现LWNM方法.LWNM的局限性在于如果壁面是运动的,LWN需要不断重新计算,而且表面固体层LWN的计算需要更多固体内部的粒子的参与,计算效率将有所降低.

4 算例验证

我们在Poiseuille流动中验证此边界条件的效果.在一个三维模拟腔中,x和y方向是周期性边界条件,在z方向布置上下两块平板,平板和流体间为无滑移边界条件,对流体施加沿着x方向的体积力.根据纳维-斯托克斯方程可以得出这个情况的精确解[24]:

其中d是两板之间距离的一半,ν是运动黏度,F是单位质量上的体积力.模拟中采用的参数设置如下: ρ=4,kBT=1.0,a=18.75,γ=4.5,σ=3.0,rc=rcw=1.0,F=0.02.模拟盒子的大小是30.0 × 30.0 × 34.0,其中上下壁面壁厚为 2.0,总共包含121500个流体粒子和16200个固体粒子.当流体充分发展后,在z方向上按0.2的厚度设置统计层,对流体粒子的x方向上的分速度进行统计,统计结果绘制于图5.

图5 LWNM边界方法统计结果和理论解的对比(Vx)Fig.5.Comparison of LWNM and theoretical results (Vx).

从图5中可见模拟结果和(16)式给出的理论解吻合得相当好,证明了LWNM能提供真正的无滑移边界条件.图6显示了在流体区域的温度和密度分布,可以看到壁面附近的区域无论温度还是密度的波动都非常小.

验证了LWNM边界方法优秀的速度无滑移控制和温度密度控制后,需要具体展示LWNM最大的优势,即对处理复杂几何边界条件时的适应能力.以下是LWNM在两个DPD应用上的使用情况.

1)具有复杂结构的超疏水表面.自然界中的很多材料常具有规律性的微结构,比如人们熟知的荷叶,表面有很多微小的突起,这些微小的突起改变了表面的亲疏水性质,使材料可以达到超疏水或者超亲水的浸润属性.这些生物材料的表面特性给了人们启发,使得在工程应用中也越来越多地使用带微结构的表面来达到特定的目的,例如除冰.当冰在固体表面凝结冰层逐渐增厚,会导致很多灾难,比如高压输电线的变重、管道的爆裂、路面的变滑,最致命的还是飞机机翼的结冰,会直接使飞机失去升力,还曾经导致了美国3407航班的空难.当冰层已经累积起来后除冰会非常困难,所以一般在冰层开始凝结时就进行干预防止凝结,常用方法有加热法和化学法,但是加热法太浪费能源,化学法又可能腐蚀表面.近年来,具有微结构的表面给除冰带来了新的研究思路,通过加工表面的微结构使表面的疏水性增强,当液态水沾到表面时不容易附着,在微小外力作用下很容易发生滚动并脱落,由此防止了在固体表面上结冰.如何通过调整材料表面的微结构而不是改变表面的化学性质来调节亲疏水性是近几年来研究的热点,研究者们设计出了各种各样的微结构,并测试它们的效果.Liu和Kim[25]还设计出了一种带二级结构的微结构表面,可以使完全浸润材料如硅,在经过表面微结构的设计后达到超疏水的表面性能,实验证明可以使任何流体在滴落后弹开而不附着.更进一步还有学者通过表面微结构的密度梯度,来控制液滴反弹的高度和方向.仅改变表面的微结构,而不用改变材料的化学性质就能带来如此巨大的表面性能的改变,这给材料工程带来了一个新的思路和研究热点.DPD等粒子方法也常被用来模拟液滴和表面碰撞的过程,由于微结构表面已经不属于简单形状的表面,传统边界条件已经不再适用,所以在以往研究中往往通过构造高密度的壁面防止流体粒子的穿透,这种方法的弊端在前文中已经阐述.采用LWNM边界方法,可以在满足各项边界要求的条件下进行模拟.图7(a)展示了文献中微结构表面的实验照片,图7(b)展示了应用本边界条件,对固体壁面粒子赋予局部法向量后的可视化结果,图中的箭头代表每个壁面粒子的法向量.

图6 LWNM边界方法统计结果和理论解的对比(温度和密度)Fig.6.Comparison of LWNM and theoretical results (temperature and density).

图7 (a)微结构表面对液滴亲疏水性的影响[26]; (b)粒子模拟中微结构表面的LWNFig.7.(a)Microstructures on surface affect the hydrophilicity; (b)LWNs (grey arrows)of surface with microstructures.

2)复杂通道.很多应用中需要用到几何形状复杂的管道,例如分叉和汇集血管中的血液流动、微流控芯片中的复杂管路等.模拟介观尺度的血液和其中的红细胞是DPD的重要应用之一,红细胞模拟的长度尺寸通常在几到几百个微米之间,非常适合介观方法DPD,此外DPD是一种粒子方法,在模拟复杂结构的流体时非常灵活,又满足守恒定律,可以从系统的状态追溯到复杂流体的相关本构方程.这些都使DPD非常适合用来进行红细胞相关的模拟研究.红细胞模拟的研究按照时间发展主要分为三个阶段: 第一阶段模拟单个红细胞,主要解决红细胞的建模问题,包括红细胞的变形和舒展,并解决单个红细胞在简单流动中的动力学行为[27]; 第二阶段模拟两个红细胞,主要关注两个红细胞之间的相互作用,包括聚集行为和解聚集行为[28]; 第三阶段也正是当前最被关注的阶段,模拟大量的红细胞在真实血管中的流动行为,比如管流或剪切流中的运动行为等[29].在前两个阶段,红细胞一般都是处于无限大不受限的流体中,或者在简单的圆管中,但是在第三个阶段就必须考虑复杂的血管形状的影响,比如血管的分叉、狭窄、以及堵塞等.对复杂的真实的血管进行建模,并施加合适的边界条件,是模拟最终能够帮助医学研究并有更广阔应用的前提之一.

图8展示了应用LWNM边界方法进行复杂血管生成的过程.首先要进行血管建模,这个模型可以是由CAD软件绘制,或是由临床扫描图像转化而来,如图8(a)和图8(b); 然后对该模型进行网格布点,所有的网格节点都将对应粒子模型中的壁面粒子,对模型提取表面层的粒子,此步可以通过对固体粒子计算固体体积分数来得到,如图8(c);最后对固体粒子计算局部法向量,形成可供DPD模拟的固体壁面模型,如图8(d).在此边界条件的支持下,血液细胞的模拟将有很多新的问题可以模拟研究.

图8 复杂血管的局部法向量的生成Fig.8.LWN generation in complex blood vessel.

微流控芯片中的流体管道通常具有很多转向和分叉汇集,我们采用LWNM边界处理方法做了一个复杂形状的管道.先通过贝塞尔曲线绘制TJU形状的路径,然后由圆形截面沿该路径形成复杂管道,接着用固体壁面粒子填充管道外的空间,再采用LWNM边界条件仅保留边界层固体粒子并计算局部壁面法向量,液滴采用MDPD模型,模型参数为 All=—40,Bll=25,rd=0.75,固-液之间的参数Asl=—12,壁面对液滴的接触角约为130°.每个液滴由1607个DPD粒子组成.模拟完成后对液滴的材质设置成水,管道壁面的材质设置成毛玻璃,再在顶部添加方形面光源,在管道下方添加背景底板,最后渲染得到图9.

图9 应用LWNM实现液滴在TJU(同济大学)形状的复杂管道中运动Fig.9.Droplets in TJU (Tongji University)shape tube with LWNM.

模拟结果表明,对这种非常不规则的管道形状,LWNM边界处理方法仍表现出了很强的适应能力,对微流控芯片中复杂管道用粒子方法进行模拟研究提供了有力的辅助.

5 结 论

本文提出一种适用于复杂几何壁面的边界条件LWNM.该方法通过对壁面粒子的周围粒子的位置进行加权计算,得到壁面局部的法向量,从而在流体粒子穿透壁面时提供可靠的位置速度修正方向; 通过定义流体粒子的周围固体体积分数来判断流体粒子是否已经穿透壁面,如果流体粒子的φ大于阈值则认为流体已经穿透固体,并根据壁面的局部法向量施加垂直壁面向外的作用力; 同时通过对固液粒子之间黏滞力的修改实现无滑移边界条件.通过Poiseuille流算例证明了LWNM边界条件可以达到速度无滑移条件,并能出色地控制壁面附近流体的密度和温度波动; 通过带微结构的表面和复杂管道的例子,展示了LWNM对多应用的广阔前景.