屈服应力流体的粗粒化分子动力学模拟研究

1、屈服应力流体的粗粒化分子动力学模拟研究摘要：非牛顿屈服应力流体广泛应用于生活和工业生产中，具有很高的研究价值。利用耗散粒子动力学(DID：Dissipative Particle Dynamics)方法来对屈服应力流体进行模拟研究%首先应用D1D方法构建了一种屈服应力流体模型，并模 拟预测了该流体模型具有Herschel-Bulkley流体的流变特性%通过数值模拟验证了模型在剪切和伯肃叶流动特性，符合屈 服应力流体的基本特点关键词：非牛顿，屈服应力流体，数值模拟,耗散粒子动力学shear and poiseuille are verified by numerical simulation,a

2、nd itAbstract:Non-Newtonian yield stress fluids are widely used in daily life and industrial production,which shows high research value.This paper uses the Dissipative Particle Dynamics (DPD) method to simulate the yield-stress fluid.A yield-stress fluid model is constructed with DPD method,and the

3、rheological properties of fluid model are consistent with those of Herschel-Bulkley (H-B).The flow characteristics of the model in accorded with the basic characteristics of yield-stress fluid.shear and poiseuille are verified by numerical simulation,and itKeywords:Non-Newtonian,yield-stress fluid,s

4、imulation,DPD(1)(1)(2)。是一个服从高斯统计分布的随机波动变量，满 足=和= (&ik&jl+&il&jk)&(tt)。C，D 和R 分别是这三种力的权重函数，这些权重函数都符合下式的形式：ri在日常生活中，很多物质在通常情况下既不是固体也不是 简单的牛顿流体，而具有非线性的流变行为，这类物质称为非牛 顿流体。其中，有一类材料只有在被施加的作用力大于某个阈值 时才会展现出流动的性能，反之则表现出固体的特性，具有这类 临界流变学特性的材料称为屈服应力(yield-stress)流体，如: 绘画颜料、泡沫材料、湿润的水泥、清洁膏、蛋黄酱和牙膏等)#-2*， 另外在石油工业中，重

5、新启动凝胶原油管道时能够准确地估计 最小压力是十分重要的。由于流动应力阈值的存在，这类材料 具有高度的非线性特点，往往对它流动的动态特性起着极其重 要的作用，具有很高的研究和应用价值。传统宏观的数值模拟方法发展起来的理论架构体系对非牛 顿屈服应力问题的预测及分析并不够准确和细致，为了能从物 质的微观给出合理的分析验证，因而需要将微观方法引入到该 问题中。这些具有独特的性质，使它们能够具有联系宏观和微观 的潜力，同时在模拟中既能兼顾计算效率和模型尺度，能够为微 流体问题的研究工作提供可参考性较高的数值指导。本文采用的是粗粒化的分子动力学模拟方法，即：耗散粒子动 力学方法，以实现构建屈服应力流体模

6、型和完成相关模拟工作。 1耗散粒子动力学方法标准的DPD方法首先由Hoogerbrugge和Koelman提 出，是MD自下而上发展而来的粒子粗粒化的方法之一。对于任 意第i个DPD粒子符合牛顿第二定律：五=!业=!dt =!，dt =其中r,、i和！分别为第i个DPD粒子的位置、速度向量和 受力之和(包括非粒子间的外部作用力)，粒子间作用力方向为 两粒子中心连线方向。另外，单个DPD粒子的质量被设为单位 1。成对的粒子相互间作用力可以分解为三个在两粒子中心连线 方向上的分力，即保守力(商)、耗散力(F：)和随机力(F：)：=$ (+)以上分力的具体形式如下所示：F= ()?=# ()( !)

7、=$ (WM(3)其中aij,#和$分别是保守力，耗散力和随机力的幅值。 /j和Vij分别是第i粒子和第j粒子的相对位移和相对速度。并且 有 / ij = r,3(r“)=从上式(4)中可知三种力都是短程形式的，即在相对距离大 于其对应的截断半径(/c)时该部分作用力归零。式(4)中有两个 可选参数，即截断半径(rc)和权重指数(n)0 一方面，针对于标准 DPD方法，每种力都应该给定相应的截断距离(rc)；另一方面， 根据Espanol和Warren提出耗散力和随机力的权重函数及 幅值应当遵循某些约束以满足涨落耗散定理。具体来说，如果其 中一个的系数和权重函数中被任意选择了，那么另外一个应该

8、 满足下式：=) * ，$ =2#kT(5)其中的kpT是系统的Boltzmann温度，在本文中统一简称 为系统的温度。为了具有高施密特数的流体动力学响应，使得随 机力与保守力采用相同的权重指数(n=1)0需要指出的是，因为 DPD方法中的随机力和耗散力满足式(5)，因此二者构成了 DPD方法内置的恒温器。本文中模拟用到的核心算法是Velocity-Verlet积分迭代算法：r(t+At) = rO + E (t) + # (At) f(t)!(t+,t)=!(t)+At !(t)!（#+#）=!（&（t+t）,!（t+t）!（t+t）=!（t）+（ #（ !（t）+!（t+t）上式中的！（t

9、）、#（t）和r,（t）分别代表着第i粒子在第t个时 刻的加速度矢量，速度矢量和位移矢量;V（t+t）代表第i粒子 在第t+ t时刻的预测速度,t取值0.01 &同时,参数是 DPD模拟方法中引入的经验参数，代表了一些随机性质的附加 效应，在本文中取值0.65。2屈服应力流体模型的构建与模拟预测为了构建具有屈服应力性质的非牛顿流体，提出了一种带 有两种类型颗粒的DPD流体模型&这里把这两种颗粒分别简称 为Pa和Pb,并且体积分数比例设置成3:7，如图1左图所示&由 于系统初始状态的粒子是规则且分成分布的，所以整体流体系 统并未达到平衡状态，即系统的温度5bT%1为了能给后期的研 究作好准备，首

10、先需要将初始状态的流体系统模拟到平衡状态， 如图1右图所示&图1屈服应力DPD流体模型初始构建分布和稳定状态的粒子图为了预测流体模型的流变性质，需要对此模型施加剪切作 用（取不同的随机力参数#为2.0、3.0、5.0）,使用Irving-Ki rk- wood模型计算得到应力张量，最终可以画出流体模型的应力- 应变曲线，如图2的上图&图2 不同#参数下（M）DPD流体模型的应力-应变图从图2中可以看到该流体存在一个明显的屈服应力So流变特性是非线性的&我们将此模拟的结果与Herschel-Bulkley流体9（以下简称H-B流体）理论模型进行对比验证， 得到了很好的吻合度，如图2的下图所示&另

11、外,H-B流体的流变特性可以用下式表达：T = T0+k%（7）式中的指数n的取值范围为0.20.8,具有剪切变稀的特 点&从图2中我们可以看到在不同随机力参数#情况下，拟合曲 线的n的取值都在以上的范围内，分别为0.31 $0.24和0.45。 3数值模拟结果 3.1剪切流验证本节中我们将对模型进行剪切作用流动的模拟验证&我们 将图3中的图像根据模拟设置分成静置形成、施加剪切破坏和 静置恢复结构等三个阶段，以便于我们做进一步的阐述分析&在 静置形成结构阶段，流体粒子经过充分地混合后达到系统平衡 状态,形成取向比较一致的拓扑结构，如图3a,在粒子Pb的体 积分数分布云图3A中也可以明显地看出，

12、各区域体积分数的 差异较大，且呈现一定的规律性&在施加剪切阶段，需要对施加 大小不同的剪切速率下的工况分别进行分析，例如图3B和b 是在施加剪切速率=0.5情况下得到的结果，可以看到各粒子 的分布结构在靠近壁面处先被剪切流动所破坏；当施加的剪切 速率增大到1.0的时候，各粒子的分布结构则被完全遭到破坏, 粒子和体积分数趋向于均匀分布，如图3C和c所示&在静置恢 复结构阶段，流体系统未受到任何的剪切作用，粒子又重新成恢 复原有的拓扑结构状态，体积分数分布亦然，这在图3D和d中 得到了充分地展示&且图4剪切作用力流场的速度剖面图：（左）屈服应力流体（右）牛顿流体图3 （M）DPD模拟中流场的体积分

13、数云图和粒子结构图的变化情况 根据剪切流场的速度剖面图，非牛顿屈服应力流体随着剪 切速率增大的同时速度剖面形状会出现一个明显的变化，具体 来说，在小于转折剪切速率的情况下，速度剖面近似地呈现一种 波浪号的特点，即在靠近中心线附近的粒子层之间没有出现剪 切流动特点，出现“块状”的流动；而在大于转折剪切速率的条件 下,其速度剖面会出现类似于图4（右）中的牛顿流体模型那样 的线性剪切流动现象&这种特殊情况的出现也可以利用屈服应 力的存在来解释，当剪切速率大于某一特定值时，流体各部分的 应力值大于屈服应力且图4剪切作用力流场的速度剖面图：（左）屈服应力流体（右）牛顿流体3.2泊肃叶流验证对屈服应力流体

14、施加多种大小不同的体积力，用以研究流 体的泊肃叶流动情况$待流体模型到达稳定状态以后，对两种类 型的流体颗粒在x轴方向上施加一个额外恒定的体积力，即 g=常数$另外，将x轴和y轴方向的边界条件设置成周期性边 界条件。经过一段时间的模拟后，可以绘制流体在x轴方向上的 速度分量随着z轴坐标的变化曲线，得到图5（上）$通过分析泊肃叶流动的速度剖面，我们可以很容易地得到速 度剖面可以分成两个部分，即抛物线形式流动部分和平台形式流 动部分$为了更加直观地展示这两种流动区域，我们绘制得到施 加的体积力为1.00（即gx=1.00）的流动模型的速度云图，如图5 （下）所示。上述的平台流动部分位于整个流域的中心线附近，且 对称分布，而抛物线形式流动部分则位于靠近壁面的区域$ 的区域大小比例也随着体积力的变化而变化。具体来说，施加的 体积力越大，