（钢铁冶金专业论文）液态纯铁的耗散粒子动力学研究.pdf

上海大学硕士学位论文 摘要 耗散粒子动力学( d p d ) 模拟方法作为国际上一种较为新型的介观模拟方法， 在模拟尺度及计算速度方面较传统微观模拟方法有较大的提高，对简单流体及复 杂流体的模拟是一种非常高效的方法，但从未见其在冶金中的应用。本文首次使 用d p d 方法对高温液态纯铁体系进行了模拟，其主要目的是初步实现d p d 方法 在冶金中的应用，并为以后实现钢液夹杂物体系介观尺度及微观、介观、宏观的 多尺度计算机模拟打下了基础。 本文在总结了液态纯铁和夹杂物研究现状的基础上，主要进行了以下几方面 的工作： 一、以分子动力学( 加) 模拟结果为基础导出了d p d 模拟参数。通过宏观 物性的对应相等，我们获得了耗散粒子动力学与现实体系时间和空间尺 度的对应关系，同时导出了耗散粒子动力学的作用参数口。 二、d p d 模拟液态纯铁体系的结果的对比。通过与m d 模拟的径向分布函数 的对比，我们发现在模拟尺度较小( 埘= 1 ，p 3 ) 时d p d 模拟可以 保留部分微观描述，而随着模拟尺度的增加，这种描述逐渐消失。 三、耗散粒子动力学模拟中模拟参数对模拟结果的影响。我们围绕d p d 模拟 内部参数卅，p 和出对模拟结果的影响进行了讨论，其合理的取值范围 为：虬 5 ，p 5 和a t o 0 4 。 通过本文的工作内容，我们基本实现了对液态铁的d p d 模拟，为后面的钢液 中氧化物夹杂的d p d 模拟研究打下了基础。 关键词：耗散粒子动力学分子动力学液态纯铁 上海大学硕士学位论文 a b s t r a c t a san e w l yg l o b a ld e v e l o p e dm e s o s c a l es i m u l a t i o nm e t h o d ，d i s s i p a t i v e p a r t i c l e d y n a m i c s ( d p d ) s h o w e dag r e a ti m p r o v m e n tc o m p a r e d 、析mt r a d i t i o n a ls i m u l a t i o n m e t h o d s ，a n dv e r ye f f i c i e n c yi ns i m u l a t i n gm e s o s c o p i ch y d r o d y n a m i c so fs i m p l ea n d c o m p l e xl i q u i d s b u tt h ea p p l i c a t i o no fd p d i nm e t a l l u r g yf i e l dh a v en o tb e e nr e p o r t e d w es i m u l a t e dh i g ht e m p e r a t u r el i q u i di r o n 、衍t hd p df o rt h ef i r s tt i m ew i t ht h ep u r p o s e o fi n t r o d u c i n gd p di n t om e t a l l u r g y t h i sc o u l dm a k eaf o u n d a t i o no fo u rf u r t h e r r e s e a r c h i n gt h a ts i m u l a t i n gi r o nm e l ta n di n c l u s i o ns y s t e mo nm e s o s c a l ea n dm u l t i s c a l e a c c o r d i n gt ot h es u m m a r i z a t i o no ft h er e s e a r c ho ft h er e s e a r c hs t a t ei nl i q u i di r o n a n di n c l u s i o n ，o u rw o r kc o u l db es u m m r i e da sf o l l o w s f i r s t l y w e g a i n e d d p dr e p u l s i o n p a r a m e t e raa c c o r d i n gt oi s o t h e r m a l c o m p r e s s i b i l i t y1 c 一，w h i c hc o u l db ed e r i v e df r o mt h es i m u l a t i o nr e s u l t so fm d t h e r e l a t i o n s h i p o fs p a t i a l - t e m p o r a ls c a l eb e t w e e nd p da n dm dh a sb e e ne s t a b l i s h e d s i m u l t a n e o u s l y s e c o n d l yw em a d eac o m p a r a t i v es t u d yb e t w e e nd p dr e s u l t sa n dm dr e s u l t s a c c o r d i n gt ot h er a d i a ld i s t r i b u t i o nf u n c t i o nw ef o u n dt h a td p dc o u l dr e s e r v es o m e d e s c r i p t i o no fm i c r o s c a l ed e t a i l ，w h i c hd i s a p p e a r e df a s f l yw i t ht h ei n c r e a s i n go f n 。 a n d p t h i r d l yw ee x a m i n e dt h ee f f e c t i o no ft h es i m u l a t i o np a r a m e t e r so ns i m u l a t i o nr e s u l t s i nd p d w ef o u n dt h a t 以 i r z e m e 川州+ 。 ( 2 4 ) 在处理大的分子聚集体系、生物大分子体系及高分子体系时，往往运用某种具体 的力场来描述体系的势能。 当动能和势能都表达清楚后，就可以建立该体系的运动方程 a h p _ _ 百 。8 h q = - 面 在笛卡尔空间进行数值积分求解。体系的温度与各原子的平均速度有关： 3 k 丁= 聊，( v ，v , ) n i = l ( 2 5 ) ( 2 6 ) 从原子的位置、速度、键接方式、各种势能函数计算出体系的总能量。然后 计算各个原子在该力场中的势能梯度，得出每个原子在分子力场中所受的力，按 照牛顿第二定律，就可以计算原子的运动行为。这是一个不断循环的过程。在全 原子模型，并没有任何限制条件下，步长约l f s ( 1 f s = l o 。1 5 s ) 。足够次数的循环将 完成体系运动方程的积分过程，而得到一个多体问题的解，和在相空间( 即原子 位置及动量坐标的空间) 中的运动轨迹。 8 上海大学硕士学位论文 2 1 2 积分方程 考虑含有个分子的系统，依照经典力学，系统中任一原子f 所受的力为势能 的梯度： e ：一v j y ：川昙+ 歹昙+ j | 导) 矿 呶ta ) io z i 因此，由牛顿运动定律可得f 原子的加速度为 e a f = 1 将牛顿方程对时间积分，预测第，个原子经过一段时间后的速度与位置。 ( 2 7 ) ( 2 8 ) 用有限差分方法求解牛顿运动方程有很多种算法，比较常用的有v e d e r t 算法、 l e a p f r o g 算法、v e l o c i t yv e r l e t 算法等【2 6 - 3 0 。 2 1 2 1v e r l e r t 算法1 2 6 , 2 7 1 v e r l e r t 算法基于体系在卜址、，时刻的位置r ( f ) 和r ( t a t ) 和在f 时刻的加速 度a ( r ) 计算体系在f + 血时刻的位置r o + a t ) 。 r o + a t ) = 2 r ( t ) 一r ( t a t ) + f 2 a o ) 其速度由式( 2 3 ) 计算 v ( f ) = r ( t + a 瓦t ) - 广r ( t - a t ) ( 2 9 ) ( 2 1 0 ) v e r l e r t 算法的优点在于对计算机内存要求小，同时缺点也较明显，由于缸2 a ( f ) 项的存在降低了计算精度，同时计算f 时刻的速度需要f + 出时刻的位置，给计算 带来不便。 9 上海大学硕士学位论文 2 1 2 2l e a p f r o g 算法0 2 8 2 9 i l e a p f r o g 算法对v e r l e t 算法进行了改进： v ( t + l a t ) = v ( 卜j 1 出) + a t a ( t ) r ( ，+ 出) = r ( f ) + 脚( h 三出) a o + a t ) ：f ( t + a t ) r n ( 2 1 1 ) ( 2 1 2 ) ( 2 1 3 ) 与v e r l e t 算法相比，l e a p f r o g 算法有包括显式速度项和计算量稍小的优点， 但它也有明显的缺陷，即位置与速度不是同步的。这意味着在位置一定时，不可 能同时计算动能对总能量的贡献。 2 1 2 3v e l o c i t yv e r l e t 算法1 3 0 l v e l o c i t yv e r l e t 算法克服了v e d e tl e a p f r o g 算法在计算体系位置和速度时的不 同步性： r ( f + 岔) = r ( f ) + ) + 丁a t 2 a ( t ) a o + a t ) ：f ( t + a t ) m ( 2 1 4 ) ( 2 1 5 ) v p + f ) = v ( f ) + j 1m a ( f ) + a ( f + 缸) 】 ( 2 1 6 ) 这种算法的优点是给出了显式速度项，并且计算量适中，目前应用比较广。 2 1 3 周期性边界条件 由于受计算机硬件的限制，即使使用巨型计算机，m d 模拟也通常只能应用 于几万至几十万个粒子的体系。但由于位于表面和内部的分子受力差别较大，将 会产生表面效应。为消除此效应，同时建造出一个准无穷大体积，使其可以更精 l o 上海大学硕士学位论文 确地代表宏观系统，必须引入周期性边界条件。把小数粒子体系看作一个中心模 拟盒子，这样就可用相对少量的粒子来研究物质的宏观物性，减少计算工作量。 模拟体系是由中心模拟盒子在各个方向上重复叠和而成，沿所有方向的模拟盒子 的影像是周期性的排列。 以二维的计算系统为例，图2 1 显示二维盒中系统的排列及粒子移动方向。图 中，位于中央的盒子表示所模拟的系统，其周围的盒子与模拟的系统具有相同的 排列及运动，称为周期性镜像( p e d o d i cm i r r o ri m a g e ) 系统。当计算系统中任一粒 子移出盒外，则必有一粒子由相反的方向移入，如图中的第2 个粒子。这样的限 制称为周期性边界条件。如此，系统中的粒子数维持恒定，即密度不变，符合实 际的要求。 为 唼善、二， 。 奄 9 。j 一 国 一囟 囝 心、黼；潼 l a 毒咚骞。b o ，一， = = j ，- 、 、一、t _ - 。一 图2 1 二维周期性粒子的移动 计算系统中分子间作用力时，采用最近镜像方法( n e a r e s tm i r r o ri m a g e ) 。如 图2 2 所示，计算分子1 与分子3 的作用力，是取与分子1 距离最近的镜像：d 盒 中的分子3 ，同样，计算分子3 与分子1 的作用时取计算系统中的分子3 与e 盒 中的分子1 。 上海大学硕士学位论文 国m 。奄 。奄 羹固” c l 分豳锄暑 t 名 l o d 暑 乃国 。 奄 鼍 o f g n - 一- 2 1 4 温度和压力控制 图2 2 粒子的最近镜像 为了使“计算机实验”迫近现实生活中的实验条件，1 9 8 0 年a n d e r s o n 首先提 出了定压、定温、定温定压三种新的平衡m d 方法【3 l 】，到1 9 8 4 年已有三种定常温 度m d 方法问世f 3 2 。4 1 。如今常用的三种控温方法是：n o s e h o o v e r 3 5 1 ，b e r e n d s e n 3 6 】， g a u s s i a n 控制【3 7 】，这些控温方法的基本原理是用一个无限大热源与所考查的体系 相耦合以维系体系的平均温度与目标温度保持一致，借用这样一个假想的物理模 型来实现对体系进行温度控制。不同的等压m d 模拟方法与等温m d 的原理是相 同的，只不过是用压力代替温度、用原子位置代替原子速度所起的作用，其具体 的方法描述可参见文献【3 6 ，3 s 。 2 1 5 物性计算 求出模拟系统中原子在相空间的轨迹后，应用统计物理原理即可得到系统的 宏观特性和结构特点。下面分别列出常用的物理参量统计形式。公式中，n 、m 、 p 、v 、1 ，、，、r 、o 、分别表示体系中的粒子数、单个粒子的质量、体系 的数密度、体系的体积、粒子速度、位置、体系温度、阿伏加德罗常数和玻尔兹 曼常数。公式中( ) 表示m d 时间平均值或系综平均。 1 2 上海大学硕士学位论文 2 1 5 1 热力学量 一般地，m d 模拟得到的首要的物性量就是热力学量，如下几个公式所示， 分别表示温度、势能、维里和压力。 2 1 5 2 径向分布函数 r = 去；刁 ，、 拈 ( 2 1 7 ) ( 2 1 8 ) ( 2 1 9 ) ( 2 2 0 ) 径向分布函数( r a d i a ld i s t r i b u t i o nf u n c t i o n ，r d f ) 或双体相关函数( p a i r c o r r e l a t i o nf u n c t i o n ，p c f ) 表征着结构的无序化程度，被广泛地应用于研究液态和 非晶态结构。径向分布函数为： 馏( r ) = 专( 若万驴+ r 一吩) ) ( 2 2 1 ) 其中g ( r ) 是，+ 办范围内找到某种原子的几率，p 是系统平均数密度，尺是 原子位置，艿是d i r a c 符号，为原子数。 2 15 3 静态结构因子 静态结构因子( s t a t i cs t r u c t u r ef a c t o r ，s s f ) 也是判断结构无序化程度的物理 量，其表达式为： “， i 、 u、 州 k 3 i 一 ， 4 = 川u y = q _ t 海大学硕士学位论文 s c 的= 专l 姜e x p c 。) i ( 2 2 2 ) 其中，n 为原子数，k 为倒格矢，5 为原子的位置矢量。对理想晶体而言， 其静态结构因子为1 ，而对于理想流体，则为0 。 2 1 5 4 均方位移和扩散系数 粒子位移平方的平均值称为均方位移( m e a ns q u a r ed i s p l a c e m e n t ，m s d ) ，其 定义式为： ，2 ( 垆万1 ( 酗n ( 旷删) ( 2 2 3 ) 其中( 0 ) 为原子，在零时刻的位移，( f ) 为原子f 在f 时刻的位移。根据爱因斯坦扩 散定律，均方位移随时间变化表征了液态金属原子的扩散行为。它与扩散系数( d ) 存在如下关系，式中c 为常数。 2 1 5 5 配位数 l i m r 2 = c + 6 d t t - - - - a o ( 2 2 4 ) 配位数( c o o r d i n a t i o nn u m b e r ，c n ) 即原子的第一最近邻原子的个数。配位 数的多少显示了该原子周围原子分布密度的大小，是最基本的结构考察因素。 2 2 耗散粒子动力学原理 2 2 1 基本原理 耗散粒子动力学( d p d ) 最早是由h o o g e r b r u g g e 和k o e l m a n 在1 9 9 2 年基于 分子动力学( m o l e c u l a rd y n a m i c s ) 和格子气自动机( l a t t i c eg a sa u t o m a t a ) 这两 种方法提出的一种新的介观模拟方法【2 2 】。 d p d 模拟方法可以看作是一种“粗粒化”的m d 方法，它省略了m d 中过多 1 4 上海大学硕士学位论文 的细节描述，使得模拟可以在更大的时间、空间尺度上进行，同时由于模拟尺度 由微观转向介观，局部的流体力学作用也成为了不可缺少的因素。在d p d 模拟中 模型粒子不再对应真实的原子或原子团，而是代表特定流体区域内质量和动量的 集合。此时粒子间的相互作用也不再由微观级别的保守力主导而引入了远程流体 力学力。同时在粒子作用力的计算方面使用了软势，即粒子所受保守力与作用距 离的关系由微观级别的多次因式转变为介观级别的线性关系，因此两粒子距离较 近时不会产生很强的排斥作用。与此同时，d p d 没有精确定义粒子的边界，因此 粒子可以相互重叠。 根据d p d 算法的定义，d p d 中粒子卜一共受三种类型的力，保守力、耗散力 和随机力。 b = 巧+ 唠+ 跨 ( 2 2 5 ) 其中巧表示保守力( c o n s e r v a t i v ef o r c e ) ，表示两粒子间的相互作用，其中嘞是f 粒 子和粒子间的作用系数，乃是f 粒子和粒子间的距离。 l 嗲= ? 1 一勺屯三芝 ( 2 2 6 ) 时为耗散力( d i s s i p a t i v ef o r c e ) ，即流体力学作用，表示两粒子间相互减慢运动速 度的粘滞阻力。其中】，为耗散因子，y 前的负号表明其作用方向与相对速度v f ，相 反。 i 垆= - y w 。( 乃) ( 屯v 扩) 屯 = ( 甄岛+ 4 ，颤) 万o r ) ，以保证 其符合g a u s s i a n 分布。e s p a f i o l 和w a r r e n 瞳3 1 认为两个权函数w d ( ) 和w r ( ) 中只 有一个可以任意选择，当这一选择固定则另一个也随之固定，两者的关系如式( 2 2 9 ) 所示： w d ( 勺) = 【w r ( ) 】2 ( 2 2 9 ) 两个常数y 和盯与温度有关：a 2 = 2 y k 丁其中t 为绝对温度，k 为玻尔兹曼常数。 当只有保守力作用时，y = 仃= 0 ，此时算法与m d 一致。 d p d 模拟过程中所有量都是无量纲的。如果( ，v ，f ) 是物理单位的长度、速 度和时间，那么在d p d 模拟中的对应量( ，1 ，力由下式给出【2 5 1 ： 厂= 厂吃， ( 2 3 0 ) 通常在使用d p d 的过程中，会对其单位进行进一步简化，通过设定m = 1 ， 乍= 1 以及基本能量单位= l ，将温度单位由r 简化为丁，同时d p d 中使用 r 作为能量单位，一般设定丁= 1 ，达到简化计算的目的。需要说明的是d p d 1 6 上海大学硕士学位论文 模型中的温度k 丁作用在于控制粒子初始速度分布( m a x w e l l 分布) 以及控制随机 力的大小。而d p d 模拟的对应的实际温度与获得作用系数和】，时的温度相等。 在后续进行的d p d 模拟中，如非特别说明，一般均使用d p d 无量纲单位。 2 2 2 积分算法 由于耗散力依赖于速度，而速度又由耗散力决定，因此在积分d p d 的运动方 程时需要非常小心。 2 2 2 1 原始模型的e u l e r 算法 最初的d p d 算法是e u l e r 算法，用以表征粒子的位置和速度 3 9 , 4 0 ： r o + f ) = r ( t ) + a t v ( t ) v ( t + a t ) - - v ( t ) + a t f ( t ) f ( t + a t ) 2 f ( r ( t + a t ) ，v p + & ) ) ( 2 3 1 ) 2 22 2g r o o t 和w a r r e n 提出的修正v e l o c i t y v e r l e t 算法1 2 4 i g r o o t 和w a 玎e n 【2 4 1 使用了修正v e l o c i t y - v e r l e t 算法，引入了参数九( 通常取 o 6 5 ) ，部分考虑了速度与力的相互依赖关系，在保证计算速度的基础上进一步提 高了精度。 i o + 缸) = 吣) + 觚) + 三( 2 c 彰 篡叫，( d ：九磐_ (232)fat + a t ) = f j ( r ( t + a t ) ，审o + f ) ) ， 、一7 v ，( h ) = v ，( f ) + 1 2 a 州，) + e ( h 2 2 3 耗散粒子动力学的发展 d p d 模拟研究目前主要有两个发展方向：一是其理论基础和算法，二是d p d 模拟的应用。对于前者研究工作主要集中于算法的优化和改进方面。最初的d p d 模拟算法采用欧拉算法，它是一个精度较差的算法。g r o o t 和n 【2 卅利用 1 7 上海大学硕士学位论文 v e l o c i t y v e r l e t 算法并提出一种改进形式进行模拟。n o v i k 和c o v e n e y 4 1 】比较了多 种差分格式以及含有非保守力体系的d p d 问题。l o w e 4 2 1 提出了一个结合分子动 力学a n d e r s e n 恒温方法的d p d 算法，该算法结构上满足细致平衡条件，此条件 表明在体系平衡时，基元转化速率一定与其逆反应相等。o t t e r 和c l a r k e l 4 3 】结合对 含有保守力的随机动力学体系近似解来改进d p d 模拟体系的恒温算法。e v a n s 4 4 1 研究了d p d 体系的传递性质如自扩散系数和剪切粘度。g r o o t l 4 5 1 引入了静电作用 使得d p d 可以用于模拟聚电解质体系。 更多的研究工作集中于d p d 模拟方法的应用。k o n g 等【4 6 】将d p d 模拟算法用 于高分子稀溶液系统，模拟了高分子回旋半径的标度性质及动力学弛豫时间等。 g r o o t 和m a d d e n 4 7 ，4 8 1 利用d p d 法模拟了由a 和b 两种链节组成的前端共聚高分 子熔体的微相分离问题，并进一步证明了流体力学对介观结构形成的重要作用。 e l l i o t t 和w i n d l e 4 9 1 利用d p d 方法研究了球形和非球形粒子的多分散混合物的聚集 状态。s h i l l c o c k 等刚研究了双分子层膜的平衡结构、横向张力及其剖面图。r o l a n d 等【5 1 】模拟了聚合物最低临界溶解温度。 还有一些工作是对非平衡流动系统进行的研究工作，如d a r i a s 等5 2 1 利用d p d 方法研究了悬浮液在一个约束系统中的流动问题。c h e n 等【5 3 】利用d p d 研究了一 个聚合物液滴稳态流动和瞬间剪切流动的性质。k i m 和p h i l l i p s 5 4 j 贝忻0 用d p d 方 法研究了流体在球体和圆柱体表面的流动问题。这是一个值得注意的动向，因为 对流动系统的模拟结果可以与基于连续力学的流体力学的计算结果进行比较，这 使得基于粗粒化的多尺度模拟方法和基于细化方法的多尺度模拟方法在流动系统 的层面上对接。 2 2 4 耗散粒子动力学的应用 d p d 方法有助于解决配方化学、高分子科学和化学工程所涉及到的复杂问题。 例如胶束形成、胶体絮状物构造、乳化、流变学、共聚物及高分子共混形态以及 通过多孔介质的流动等研究。它可以模拟真实实验处理条件( 压力、温度、处理 时间等) 下聚合物或胶体溶液的化学形态、微观形貌、相分离以及流变性等，为 聚合物复杂流体的研究提供了强大的理论预测工具。 1 8 上海大学硕士学位论文 与要求详细描述原子的分子模拟相比，d p d 使用比原子尺度大的基本单元， 通过动力学模拟来确定这些模型的结构、性质和动力学演变过程。与分子模拟相 比，d p d 在大得多的空间和时间尺度上描述介于微观和宏观之间的介观体系。因 此d p d 方法的应用吸引了更多人的注意。 目前d p d 已经成功应用于嵌段共聚物、表面活性剂、胶体等体系的研究。k o n g 等酬将d p d 模拟算法用于高分子稀溶液体系，模拟了高分子回旋半径的标度性质 以及动力学弛豫时间等。g r o o t 和m a d d e n 【5 叼利用d p d 方法模拟了由a 和b 两种 链节组成的嵌段共聚高分子熔体的微相分离问题。发现当聚合物的对称性减小时， 平衡的介观结构按照层状、孔形层状、六边形柱体和胶束次序变化。平均场理论 预测的稳定螺旋相结构在模拟中发现朝六边形孔形层状相转化。他们认为d p d 中 的流体力学作用对介观结构的演化起着重要作用。g r o o t 和m a d d e n 等【5 7 】结合布朗 动力学进一步研究了嵌段共聚物微相分离来证明流体力学作用所起的重要作用。 他们发现d p d 模捌舅荡彪五殆撰拟缪绣的口，记，- 妨露巧行为。g r o o t 5 8 1 利用d p d 方法研究了聚合物表面活性剂在水溶液中的聚集情况，并利用w i d o m 插入粒子法 研究了表面活性剂的化学势。y a m a m o t o 等【5 9 ，删利用d p d 模拟方法研究了水溶液 中两亲分子自发形成囊泡的现象。s h i l l c o c k 等【5 川研究了双分子层膜的平衡结构、 横向张力及其剖面图。r e k v i g 等【6 1 】研究了单尾和双尾表面活性剂在油水表面的情 况。 同时国内研究者在d p d 方面也开展了卓有成效的应用性工作。刘扬等蚓用 d p d 方法研究了聚合物材料的玻璃化转变温度性质以及界面张力对共混物形态结 构的影响。苑世领等研究了聚合物与表面活性剂的相互作用6 3 1 以及水溶液中聚合 物聚集的动力学变化过程【删。魏克成等则应用d p d 方法研究了嵌段共聚物胶束 体系的微相分离性质。 在商业应用方面，d p d 模拟方法已有商业软件可以购置，a c c e l r y s 公司的m s m o d e l i n g 软件已经包含了d p d 模拟模块。但限于d p d 技术本身还不是非常成熟， 而且应用领域也还比较有限，使得相应的商业软件中的模块的开发应用还远未达 到成熟。 1 9 上海大学硕士学位论文 2 2 5 耗散粒子动力学优势及局限 相比传统的微观模拟方法( 如) ，d p d 的优势在于粒子粗化和采用软势使 得模拟体系的时间尺度和空间尺度都有较大提高，且在精确性方面达到了较好的 水平，在计算效率和计算费用间做到了较好的平衡。 与此同时其劣势也较明显，由于粒子粗化使得d p d 模型放弃了对原子尺度的 细节描述，因此无法模拟体系的微观结构方面的性质。 2 2 6 耗散粒子动力学模拟步骤 到目前为止，d p d 技术的一般模拟步骤包括如下内容： 1 确定模拟体系； 2 设计d p d 模型，将真实模拟体系映射至d p d 模型，确定在d p d 模拟中粒 子的结构，对模拟体系进行d p d 单位归一化处理； 3 r 确定a o 、丫、a 、m 、席模拟参数； 4 确定合适的时间尺度和空间尺度； 5 将相关参数输入模拟软件开始模拟； 6 模拟结束得到模拟结果； 7 分析模拟结果，提取有用信息，将结果单位由d p d 单位转换至真实模拟体 系单位。 2 2 7 通过实验可测性质拟合d p d 中相互作用参数 由于d p d 使用了一套自有的单位体系并使用了“软势 作为其作用势，使得 d p d 体系与真实体系的对应以及作用参数a 的获取成为了d p d 模拟中的核心问 题。 2 2 7 1g r o o t 和w a r r e n 法1 2 4 i 在d p d 模拟中最核心的问题是模拟体系的作用力与具体研究体系性质间的映 2 0 上海大学硕士学位论文 射关系。在这类研究中，g r o o t 和w a r r e n 2 q 的工作目前来说是最具指导作用的。 g r o o t 和w a r r e n 对具有形如式( 2 2 5 ) 的作用力的体系进行了模拟。作用参数a 分别选取1 5 、2 5 和3 0 ，模拟数量密度p = 1 1 0 。将( p p k b t ) ( a p 2 ) 与p 的关系 作图，发现当密度较高时( p 2 ) 所有体系落在同一曲线上，且该比值接近一个 常数a = o 1 0 1 + 0 0 0 1 ( 见图2 3 ) 。 p 图2 3 ( p - p k n t ) ( a p 2 ) p 曲线口4 1 因此在较高密度( p 2 ) 下d p d 体系的状态方程可表示为： p = p k b t + a a p 2 ( 0 【= 0 1 0 1 + 0 0 0 1 ) 上式对应的等温压缩率为： k 一= 上k b t ( 望8 n 1 = 1 + 2 , m p 丁 ( 2 3 3 ) ( 2 3 4 ) 通过式( 2 3 4 ) 我们获得了体系压缩率与d p d 作用系数及数密度的关系，因此我们 只要获得了模拟体系的压缩率，并设定了数密度和r 便可算出a 。下面以水为例 子做一个简单介绍： 2 l 上海大学硕士学位论文 有： 在一个大气压和室温( 3 0 0 k ) 条件下，水的压缩率约为1 6 。由式( 2 3 4 ) 对水 盟7 5 k 丁 ( 2 3 5 ) 如果选择模拟体系的质量密度为p = 3 ，k r = 1 ，则水的作用参数为口一，= 2 5 k 丁。 虽然g r o o t 和w a r r e n 给出了计算a 的方法，但是并没有给出d p d 体系与真实 体系的对应关系，若想通过式( 2 3 0 ) 来对应我们还需设定d p d 粒子的质量及截断 半径，因此其实用性并不高。 ( ；r o o t 和w a r r e n 进一步将上面的d p d 模型自由能表达式与高分子体系的 f l o r y h u g g i n s 理论进行比较，得出不同高分子链节间的力参数口 ，和f l o r y h u g g i n s 参数勋间有如下关系： n q2a i i + 3 2 7 x g ( 2 3 6 ) 从式( 2 3 6 ) 可知，只要知道了体系的f l o r y - h u g g i n s 参数即可获得d p d 模拟的 力参数，进而用d p d 进行介观结构的模拟，这就为d p d 模型与实际体系间提供 了一种映射关系。对于一个具体的体系，可以有多种方法获得f l o r y h u g g i n s 参数。 最直接的方法是用f l o r y h u g g i n s 模型关联高分子体系的相平衡或其它热力学性质 的实验结果从而得到f l o r y h u g g i n s 参数】，也可以用分子力学的方法获得粗粒化 粒子体系的混合能进而获得f l o r y h u g g i n s 参数【6 3 ，d 7 1 ，还可以利用分子力学结合 m o n t ec a r l o 模拟来获得f l o r y h u g g i n s 参数【6 引。还有一种称之为i n v e r s em o n t e c a r l o 的模拟方法也可以用来粗粒化实际分子而用于d p d 模拟工作1 6 9 】。 g - r o o t 和w a r r e n 的这套参数推导方法优点在于计算简单，f l o r y - h u g g i n s 参数 可通过多种方法获得，而得到f l o r y - h u g g i n s 参数就可直接算出保守力参数。但是 其弱点也非常明显。首先由于d p d 软粒子模型和f l o r y h u g g i n s 格子模型有着很 大差别，而且软粒子势能的描述既可能是空间上的粗粒化也可能是时间上的平均， d p d 作用参数和f l o r y h u g g i n s 参数的这种联系并不是非常严格的，而仅仅是根据 上海大学硕士学位论文 两个体系的某种宏观性质的类比得到的；其次在他们模拟的体系中粒子间的势能 为纯排斥势能，与实际情况也是有一定偏差的；第三，f l o r y h u g g i n s 理论是针对 常温高分子溶液而设计，并不适用于我们要模拟的高温钢液夹杂物体系。 2 2 7 2k e a v e n y 方法1 7 0 l 此法以g w 法为基础，通过将d p d 体系中质量浓度、粘度和压缩率对应相 等而得到d p d 模拟参数a 与体系压缩率，以及d p d 与m d 之间时间和空间尺度 的对应关系。我们将d p d 粒子和m d 分子的关系定义为：m = 虬聊，式中m 为 d p d 粒子质量，m 为m d 中分子质量，卅为粗化系数，即一个d p d 粒子对应的 m d 分子数量。 通过将d p d 和im d 体系中的质量密度m p d p d 和聊对应相等( p o t , d 和p 加分 别为d p d 和m d 体系的数密度) ，得到长度单位的对应关系： ：出血) ；。 p ( 2 3 7 ) 式中为d p d 截断半径，o 为m d 长度单位，麻d 数值与d p d 数密度( 细d ) 相同但单位为a - 3 ，为m d 数密度。 通过动力粘度对应相等，获得时间尺度的对应关系如式( 2 3 8 ) ：vdpddpd仁) z t 2 ( 三) 。t v m b 6 ( 2 3 8 ) 通过将等温压缩率k - 1 对应相等，获得作用系数a 。等温压缩率定义式如下： ，c = 寿t ( 刍o n r 、 “ 将式( 2 3 3 ) 带入上式，可以得出对于d p d 系统的等温压缩率为： ( 2 3 9 ) l 海大学硕士学位论文 k 一= 古c 去，r 2 忐c 协等， 通过上式我们只要确定k 一，卅和肪阳便可获得作用参数a 。 其中1 c - 1 可通过k i d 模拟获得，在l v i d 中压力计算式为： 脚”古 5 时，出 现在o 8 咫的峰逐渐分化为两个峰( 埘变化时，0 8 吃对应的实际长度也在变化， 对应关系见表4 1 ，且埘越大趋势越明显，正如上文所说，由于卅越大时相互作 用系数口越大，也就是说在模拟体系大小不变的情况下，虽然模拟数密度并未变 化，但是过大的口使得模拟体系显得“拥挤”。如此大的作用力下原有的粒子排列 已无法维持，最终形成了新的排列结构。这个过程反应在径向分布函数图上就是 较小虬时的单峰随帆变大逐渐分峰的过程。k e a v e n y 的研究认为d p d 模拟中心 取值上限为5 ，此时其结果m d 模拟结果较好符合【_ 7 0 l ，我们的结果验证了其观点。 3 7 上海大学硕士学位论文 i n t e r b e a dd i s t a n c e ( d p du n i t s ) 图4 6 埘= 1 ，5 ，1 0 ，1 5 ，2 0 时的径向分布函数 4 2 2 2 a r m 对扩散系数的影响 不同虬时的扩散系数如图4 7 ( a ) 所示，扩散系数随心变大而减小，为精确 s i m u l a t i o nt i m e ( d p du n i t s ) 图4 7 ( a ) 虬= 1 ，5 ，1 0 ，1 5 ，2 0 时的扩散系数 3 8 兰c：o止。一参|f_竖o 上海大学硕士学位论文 研究其变化程度，我们选择【啪= 1 9 0 时的扩散系数对朋做图，其关系如图4 7 ( b ) 所示： 卫 芒 3 o 正 。 蚤 ； 历 ： 1 1 = a 图4 7 ( b ) 扩散系数随虬变化曲线 从图中可以看出虬的增大导致d p d 粒子的尺寸变大，此时模拟体系内部的“拥 挤”影响了粒子的运动能力，大部分动能消耗于耗散力，从而导致扩散系数的降 低。同时由扩散系数定义式( 4 2 ) 可知，扩散系数与蠢成反比，而式( 2 3 7 ) 定义了肼 与露成正比，因此我们的结果是合理的。 4 2 3p 对径向分布函数的影响 由于p 的改变会影响到与现实体系的对应关系，因此我们将径向分布函数的 横坐标进行了处理。径向分布函数的变化情况见图4 8 ( 虬= 1 ) 和图4 9 ( 虬= 5 ) 。 由于式( 2 3 3 ) 的关系，为保证精确度我们选择了p = 2 ，3 ，5 三个值作为研究的对象。 不同p 对应的a 及其与m d 体系的对应关系见表4 2 。 3 9 上海大学硕士学位论文 图4 8 密度对径向分布函数的影响 从图4 8 中我们发现随着p 的增大，径向分布函数的峰位逐渐后移，峰型逐渐 平滑，当p = 5 时峰基本消失。可以看出在坛= 1 时，p = 2 和p = 3 径向分布函数峰位 非常接近。结合表4 2 分析，此时d p d 体系无论在时间尺