umeriealSimulationofPolymerFlowwithDissiPativePartieleDynamiesMethodDissertationSubmittedtoZhengzhouUniversityinPartialfulfillmentoftherequirementforthedegreeofMasterofEngineeringbyZhangWenchong(MaterialsProcessingEngineering)DissertationSuPervisor:ProfessorChang一yuShenProfessor矶/e1CaoMay，2009 N 郑州大学硕士学位论文摘要聚合物流体在宏观尺度上观测到的物理现象，是由流体的微观结构决定的，而有些流体特征依赖于高分子的聚集态结构而不是最细微的分子结构，这就要求采用的研究方法能达到介观尺度。虽然元胞自动机和格子一波尔兹曼方法都属于介观尺度的数值模拟技术，但是它们比较难于处理复杂流体;分子动力学由于所需计算量太大而难以实现。耗散粒子动力学(DPD)是一种新兴的介观尺度数值模拟技术。在DPD系统中，其基本的单元是一些被称为粒子的动量载体，这些粒子代表的是大量分子的集合体。DPD方法的优越性在于能够容易地模拟流体系统，提供流体的分子取向、速度等微观信息。本文利用耗散粒子动力学这种方法对流体的动力学行为进行研究讨论，研究内容主要包括:1.系统讨论了DPD方法的基本理论。包括其基本思想、详细算法以及各种参数的选取等。2.应用本文方法模拟了简单流体泊肃叶流动和方形通道内流体的三维充分发展流动，并与理论分析解、经典文献结果、Flueni软件的模拟结果进行对比，验证了算法的正确性。3.对聚合物流体进行基于DPD方法的模拟，在介观尺度上计算了聚合物流体流变性能参数;结合经典流体力学方程和DPD理论方法介绍了聚合物流体的模拟过程的宏观-微观祸合分析方法。本研究工作采用了耗散粒子动力学模拟了流体流动行为，有利于更加深刻地了解流体的流变形能。关键词:耗散粒子动力学;聚合物流体;宏观一微观祸合;聚合物流变性能郑州大学硕士学位论文ABSTRACTThebehaviourofPolymerliquidobservedinmaero一seale15deeidedbyitsmiero一strueture.SomePolymerliquid’5charactersdePondonitsmoleeuleassemblingstructureratherthanmoleculestrueture.ThisrequiresthesimulatingteehniqueshavetobemesoseoPieseale.CellularautomataandBoltzmannmethodarethemesosealenumeriealsimulationteehniques，buttheyarehardtohandleeomPlexfluids.Moleeulardynamiesmethodinvolvestoomuehcaleulation，itdoesnotworkeither.DissiPativePartieledynamies(DPD)，whosebasieunit15knownasthemomentumveetorPartielesrePresentingalargenumberofelementsofeolleetivebehavior，15anemergingmesosealenumeriealsimulationteehnique.ThesuPerioritiesofDPDmethodarethatiteaneasilysimulatethefluidsystemandgetmoremeso一sealeinformation.InthisPaPer，thedynamiebehaviorofPolymertobediseussedwithdissiPativeParticledynamiesmethod.TheconduetofPaPerinelude:1.DiseussthesystematieProeessthathowtouseDPD.Ineludingitsbasieidea，adetailedalgorithlnandthemethodthathowtoseleetvariousParametersand50on.2.SimulatethesimPlePoiseuillefluidflowwithDPDmethod.ComParethesimulateresultswiththetheoretiealanalysisofsolutionandresearehresultsbyotherstoverifythecorreetnessoftheProeedure.Simulateresultofthefluidflowinsquaremierochanneleoineidewithsoftwaresimulationresult.3.SimulatethePolymerfluidflowwithDPDmethod.CaleulatetherheologiealProPertiesofPolymerfluidinmesoseoPicseale.CombinatetheelassietheoryoffluiddynamiesequationsandDPDmethodtointrodueethemaero一mieroeouPledsimulation.Inthisstudy，usethedissiPativePartieledynamicsmethodinthesimulationofPolymerfluidflowbehavior，ProvideamoreProfoundunderstandingoffluidflowrheologiealPr0Perties.Keywords:dissiPativePartieledynamies:Polymerfluid;miero一maeroeouPled;PolymerflowtheologiealProPerties 郑州大学硕士学位论文目录摘要....................……，........................................................................................................……ABSTRACT..............................................................................................……，....................……11目录...................................................................................................................................……第一章绪论1计算机模拟技术....................................................................................................……12数值模拟在聚合物材料加工方面应用................................................................……23数值模拟方法的比较及选题意义....................……，.............................................……21.3.1传统方法......................................................................................................……31.3.2蒙特卡罗方法.........................……，..............................................................……31.3.3分子动力学方法......................................................……，.............................……41.3.4格子一波尔兹曼方法.....................................................................................……61.3.5选题意义......................................................................................................……74耗散粒子动力学国内外研究现状........................................................................……75课题主要内容.......................................................................................................……，96本文结构................................................................................................................……9第二章耗散粒子动力学理论基础...............……，.....……，..…，.，，.，…，.，.，.，.……，...........……，，.…102.1DPD基本理论......................................................................................................……102.1.1标准DPD方法..........................................................................................……102.1.2改进DPD方法........................................................................……、...........……142.2模拟基本原理......................................................................................................……152.2.1无量纲方法................................................................................................……152.2.2聚合物FENE型结构.................................................................................……172.2.3粒子初始位置............................................................................................……192.2.4粒子初始速度............................................................................................……202.2.5边界条件及处理方法.......................................……、..................................……21 郑州大学硕士学位论文2，2.6粒子相互作用计算方法............................................................................……242.3DPD模拟过程.................................................................................................……，.…252.3.1给定初始条件.............……，...............................................……，..................……252.3.2重要参数选择................................……，.........................……，.....................……252.3.3系统的运行.，................................................……，.......................................……272.3.4宏观物理量的计算....................................................................................……292.3.5模拟程序.................................................................................................……、二302.4本章小结..............................................................................................................……34人含—二亡己功线涅，户去才去‘户古二卜」褪雪七、!，尸牙弓二二早l司井三切乙件训L不刃魂矢j扒..............................................……，..........……，......................……JJ3.1泊肃叶(Poiseuille)流模拟..............................................................................……353.1.1泊肃叶(Poiseuille)流....................................................……，.................……353.1.2模拟方案....................................................................................................……363.1.3数据处理....................................................................................................……383.1.4结果讨论....................................................................................................……413.2方形通道内流动的模拟......................................................................................……433.2.1模拟方案....................................................................................................……433.2.2结果讨论....................................................................................................……433.3本章小结..................................................................................................……，.....……44第四章聚合物流体流动模拟............................................................................……，.....……454.1剪切流动模拟............……，...................................................................................……454.1.1模拟方案....................................................................................................……454.1.2模拟结果及数据处理.........................................................……，................……4.2FENE链构象变化................................................................................................……484.3聚合物流变性能...........................................................……，................................……504.3.1聚合物的粘性..............……，..........................................……，......................……504.3.2聚合物的弹性............................................................................................……514.4宏观一微观尺度的祸合.........................................................................................……524.4.1流体力学基本方程....................................................................................……52 郑州大学硕士学位论文44，2宏观一微观祸合分析...................................................................................……534.4.3祸合结果讨论.........................................................……，............................……554.5本章小结..............................................................……，.............……，.....................……55竺石吝匕壮匕鼠胡《勺二川工‘仔‘，山、二口一勺夕p屯曰巳。。。。。。。。···········……。。·。。。。。。·。。。。·。。·。。。.…。。…。.…，。。。。“二，二，··.……，“··.……，··.……，.……J，参考文献............................................................................................................................……58致谢................................................................................................................................……60郑州大学硕士学位论文第一章计算机模拟技术绪论随着计算机硬件、软件行业的迅猛发展，计算机模拟作为一门新兴的交叉学科，有着越来越广泛的应用价值。计算机模拟可节省人力、物力，缩短实验周期，先验实验的可行性，预测实验的结果。还可实现实验条件下达不到的极限条件的模拟，如高温等离子体、核反应堆模拟等。如今，计算机模拟己应用到几乎所有的应用科学领域，如材料科学、生命科学、天体学、气象学以及建筑、工程学等方面。计算机模拟方法与理论和实验一样，都是认识自然的主要途径和工具。理论分析通常总是建立在高度简化的假设模型之上。而对于复杂的体系，实验有时候往往很难控制。计算机模拟可弥补这方面的不足，可作为它们之间的桥梁。在材料科学中最活跃的一个领域就是制备、表征和应用高分子材料。由于高分子的化学复杂性、组成多样性，能与其它有机及无机材料相互混和，因此高分子聚合物有广泛的材料应用前景，在国计民生中具有越来越重要的地位。高分子材料的日益增长和广泛应用是惊人的，可以说没有高分子材料的现代生活是难以想象的。世界上每一个工业化的国家都在积极开展这方面的研究工作，这使得高分子科学得到加速发展。在工程技术方面，工业高分子合成上的一系列成就开创了国际化学工业的一个新的分支，它专门致力于生产及应用高分子材料，特别是塑料、橡胶、纤维、涂料和粘合剂等。聚合物材料最主要的应用价值在于它优良的力学性能，这是由高分子内在的长链结构所决定的。通过改变诸如高分子的化学组成及其形态，我们可以得到大量的性能各异的新型高分子材料。由于高分子是复杂的体系，高分子科学本身正在不断的前进中，尚有许多需要解决、发展和突破的研究课题。除了实验手段的不断发展完善外，关于高分子链性质的非平衡态统计理论和标度理论以及近代物理学中诸如自洽场、重整化群方法、量纲分析等[l一81也正应用于高分子体系的问题处理中。近年来，计算机在高分子研究中的应用发展很快，建立在严格的量子力学及统计力学基础上的计算机模拟方法，不但使难以解析的许多理论处理中的数学问题实现数值计算，还可以直接模拟高分子科学领域中的某些物理过程和化学过程。郑州大学硕士学位论文三面面面面面面面面面面百面面蔺面从实用的角度来讲，要获得具有特殊性能的高分子材料，常需要进行大量的实验，然而这需要大量的资金和时间，并且取决于实验技术人员的经验。与实验科学相比，计算机模拟方法有着不可比拟的优越性，是帮助我们理解高分子的结构性质、预测材料的性质、行为的理想选择。例如，它可以准确预测通过改变高分子的立构规整度而导致的材料性能的变化，而这在实验中由于条件是难以实现的。.2数值模拟在聚合物材料加工方面应用聚合物成型加工过程的计算机模拟主要是利用高分子材料学、流变学、传热学、计算力学和计算机图形学等基本理论，建立塑料成型过程的物理和数学模型，构造有效的数值计算方法，实现成型过程的动态仿真分析，使人们对塑料成型过程的认识从宏观进人微观、从定性进人定量、从静态进人动态，为优化模具设计和控制制品成型过程以获得理想的最终产品提供科学依据和设计分析手段。在注塑成型中，对填充过程进行数值模拟可以预测实际注射过程中可能出现的缺陷，尽快优化模具结构设计、调整工艺参数和有针对地性制订技术解决方案。以上所述是宏观尺度的模拟。众所周知，材料内部的微观与介观结构对于制品的性能起了决定性的作用。微观结构主要取决于材料的选择，由于其松弛时间很短，加工成型过程对其影响不大(在不考虑化学反应时)。介观结构(如聚合物熔体大分子链构象)尺度较大，松弛时间较长，在流场作用下处于非平衡状态，在冷却时易被冻结下来，同时引发内应力，影响结晶过程。因此，如何描写聚合物大分子链形态在加工过程中的变化是当前聚合物加工过程数值模拟的新的研究热点。.3数值模拟方法的比较及选题意义在我们的研究工作当中，很多实验对实验器材的精度要求很苛刻，细微的误差就会导致结果的失真，同时理论分析方法的应用范围也非常有限，很多复杂情况即使采取一些假设，也不能顺利地解决问题。这时，计算机模拟方法就有了一定的应用空间。它不仅将模型与理论预测联系起来，也将模型与实验结果联系起来，因此，人们也常把计算机模拟称为“数值实验”191。为了研究高分子在各种尺度下的性质，在微观和宏观方面，人们均已经建立起一些郑州大学硕士学位论文较成熟的模拟方法，微观方面的模拟方法，如:分子动力学(MoleeularDynamies，MD)、MonteCarlo(MC)、布I胡动力学(BrownianDynamies，BD)等。宏观方面的模拟方法，如:基于流体力学的有限元方法(Finite一ElementMethodbasedonF一uidMeehanies，FEMBFM)。但由于高分子的长链结构，使其物理性质在时间和空间尺度上都表现为多尺度性，人们时常希望不仅取得高分子在宏观方面的性质，而且还能获得其微观方面的重要信息，因此，介观方面模拟方法应运而生，如:耗散粒子动力学(Dissipativepartieleoynamies，oPo)和格子一玻尔兹曼方法(LattieeBoltzmannMethod，LBM)，动态密度泛函理论(DynamieDensityFunetionalTheo以DDFT)以及场论聚合物模拟(Field一theoretiepolymersimulation，FTPS)等。介观模型方法的特点就是在快速分子尺度和慢速宏观尺度之间架起联系的桥梁。在目前除了一些功能强大、技术成熟的大型软件己用于科学研究、工程技术上的常规计算问题，由于发展时间相对较长，我们将其归类称为传统方法外，最常用的几种模拟方法主要是分子动力学，蒙特卡洛和格子一玻尔兹曼方法。这些方法分别在不同的领域表现出自身的优势，下面逐一将各种方法与DPD方法进行分析和对比，从理论基础、应用领域及计算效率等方面分析各自的优劣，同时明确了耗散粒子动力学方法在数值模拟方法中的地位。1.3.1传统方法传统方法主要是应用软件，基于N一S方程的连续介质流动模拟方法，如有限元方法、有限容积法和有限差分法等。传统方法技术相对成熟，应用范围很广，在常规尺度下可以充分表现出优势，但对处理像涉及到微观结构的多相流问题，传统方法在界面处理上常常离不开人为的干预。1.3.2蒙特卡罗方法蒙特卡罗方法，又称随机抽样或统计试验方法，是以统计力学为基础的一种数学实验方法，将这种方法应用于数值模拟，就可以方便地解决很多复杂问题。美国科学家Metropofis首次对复杂的化学体系实现蒙特卡罗数值模拟[’“]，开创了现代Montecarto模拟的先河。蒙特卡罗方法的;基本思想是:当所要求解的问题是某种事件出现的概率，或者是某 郑州大学硕士学位论文个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。蒙特卡罗方法是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以归结为三个主要步骤:构造或描述概率过程;实现从己知概率分布抽样;建立各种估计上巨。蒙特卡罗方法的主要特点:直接追踪粒子，物理思路清晰，易于理解;采用随机抽样的方法，较真切的模拟粒子输运的过程，反映了统计涨落的规律;最大的优势在于它不受系统、多因素等复杂性的，是解决复杂系统粒子输运问题的好方法。例如，现实世界中许多问题的影响因素很多，问题的解决都归结为问题。一般的数值分析方法在处理这类问题会因求解所需的计算量随着问题的维数的增加呈指数状增加，而蒙特卡罗方法对维数不敏感，可以很有效率的实现模拟。因此，MC方法不足之处在于解的收敛速度不理想，而且概率误差对随机样本的要求很高。1.3.3分子动力学方法 不同于MonieCarfo这种基于概率的随机性模拟，分子动力学模拟是另一种处理复杂体系的计算机模拟方法[9，”一，41。MD方法中心思想就是计算构成系统的分子或原子之间的相互作用力，然后采用经典的牛顿运动方程求解，对于量子效应不明显的体系，可以得到比较满意的结果。分子动力学方法对于所研究的体系，考虑其中每一个分子或原子间的相互作用势，可写为: 。·艺vl(。)+艺艺vZ认，r，)+艺艺艺v3认，rj，。)…J>lk>j(1.1)然后将这个体系采用Lagrange形式的运动方程描述为:一dt}。竺}}由:)又Uqi/、r二)‘耐.—}一}—}二U、‘’叼(1.3) L二K一V郑州大学硕士学位论文‘=万六(+式x/，·/)zq，共扼的的定义为: (l，4)L称为Lagrange算子，q和p分别为广义坐标和广义动量，广义动量p，与广义坐标况Pi=一1一日q*在直角坐标系中，上式可化为:(1.5) m，‘=关，关二v。L=一v。V(1.6)如果采用与Lagrange形式的运动方程对应的Hamilton形式的运动方程描述即:日万予;|厂L.|qi二二厂一印， (1.7)其中H称为Hamilton算直角坐标下，Hamilton运动方程可化为:‘Il|t| 二mi科一V二拼P关(1.8)不管采用哪一种形式表达的运动方程，都将分子间势能、分子间作用力以及分子的位置、速度、加速度联系起来。这样，求解体系随时间的演化的表现，只要求解3N个二阶微分方程(Lagrange形式)，或是与之等价的6N个一阶微分方程(Hamilton形式)即可。与MC方法一样，MD方法也是要产生一个描述体系状态的统计系综的抽样。如果感兴趣的仅仅是静态的平衡性质，就可以生成一个由平衡态构成的系综。但若要研究动力学性质或非平衡态的性质，那么就要使用动力学方法来生成相空间的轨迹。不管是哪一种状态的系统，要想得到正确的对体系的描述，就必须保证得到精确的微观分子模型以及相应的高精度的原子间相互作用势。同时，也要做到对构型空间的充分取样。郑州大学硕几l:学位论文面面面百面面面函面面面函面面面面分子动力学方法可以获得研究领域所感兴趣的微观特性，如原子、分子间相互作用等，但是这种优越性由于目前计算机速度的而仅限于纳米纳秒量级范围内，在毫秒微米量级即介观尺度上的应用还指望将来计算机速度的提高。例如高聚物体系，生物大分子，胶体等“复杂流体”，动力学过程相对缓慢，即使用现有的最快的超级计算机，对包含百万个原子的体系进行分子动力学模拟，往往要几个月的时间，这正是MD方法的弱点，计算效率很低，耗费时间巨大。1.3.4格子一波尔兹曼方法格子一波尔兹曼方法是以流体力学为基础发展起来的一种介观尺度模拟方法，由格子气自动机(LatticeGa:Automata，LGA)发展而来，采用离散粒子在一定的格子上运动和相互作用实现对流体特性的模拟。格子气自动机LGA模型[l51与传统的把流体作为连续的介质来描述的数值计算方法截然不同。它将流体及其存在的时间、空间完全离散，把流体想象成由大量只有质量没有体积的流体粒子组成，所有这些粒子同步地随着离散的整时间步长依据给定的简单碰撞规则在网格节点上相互作用碰撞，沿网格线在节点之间运动，碰撞规则遵循质量、动量和能量守恒定律。它认为流体运动的宏观特征是由微观流体粒子相互碰撞而在整体上表现出来的统计规律。它既可以描述复杂边界流动问题的宏观流动特征，又能得出流动的微观细节。最初将流体场划分成正方形网格的HPP二维模型，单位质量的粒子在正方形的格子线上移动。由于正四边形的各向异性，不能正确地推导出Navier-Stokes方程，无法充分反映流体的特征。此后提出了二维正六边形FHP模型，FHP模型宏观上的行为符合Navier一stokes方程。后来由于模拟相变的需要，在FHP模型基础上，建立以一种流体的气液两相的粒子为研究对象，相邻粒子的碰撞规律与FHP模型完全一致，另外增加了相距数格子粒子之间的相互作用引力、斥力的LGliquld一gas模型，并设计相应的碰撞规律来实现液体相变过程，使格子自动机有了较大的发展并趋于完善。由于格子气方法的节点用一系列的0或者1来表征每个速度方向上是否存在粒子，所以对于格子气的运算只涉及到一些布尔运算，然而这样存在不少的固有缺点，如粒子平衡分布为Fermi分布，从而导致非伽利略不变性，既N一S方程对流项前面有一个依赖于密度的系数为1，压力依赖于宏观速度，随即统计噪声比较大。 郑州大学硕士学位论文为了克服上述缺陷，由Boltzlnann方程去替换格子气自动机，简单的用分布函数去代替布尔变量n，求f的方法是根据碰撞规则表，这样就把原来的整数运算变成了实数运算，解决了随机噪声问题。这就是我们前面提到的格子玻尔兹曼方法[’“一”]。但是LB方法的不足就是人为的设置格子，在处理复杂边界时不够灵活;另外了粒子之间的自由运动，这就需要更多的处理方法来使其更加合理。1.3.5选题意义到此为止，我们介绍了目前数值模拟中比较常见的模拟方法，其中格子玻尔兹曼方法和后面将详细介绍的耗散粒子动力学被称为介观尺度模拟方法，目前在解决空间、时间上微米、毫秒尺度系统的流体流动问题方面存在着极大的优势。DPD和MD方法都是无网格方法，通过求解粒子的运动轨迹确定流体系统的运动规律，但由于每个DPD粒子包含了大量分子或原子，其运动的特征时间和特征长度比MD粒子大得多，因此DPD方法可以处理较大时空尺度问题，并且与分子动力学相比，耗散粒子动力学所需要的计算量较小。DPD方法比LBM方法更灵活，流体粒子的运动方向不再局限于格子网格线。另外，DPD方法还可以处理宏观连续介质模型如Navier一stokes方程不再适合的微流动、复杂流动等粒子特性明显及壁面有滑移的情况。 而且，一般对聚合物溶液模拟采用蛇行管模型，即Lennard一JonesModel，但这个模型是基于电子层面的，计算量非常大，而且各个蛇行管不能相互穿越，DPD模型可以克服这两个问题。耗散粒子动力学具有分子动力学灵活自如以及格子一波尔兹曼方法处理介观尺度问题的特点，因此本课题拟将DPD法应用于聚合物流体流动数值模拟方面，即基于DPD模型进行聚合物流体流动数值模拟。.4耗散粒子动力学国内外研究现状 Hoogerbrugge等人l’8]于1992年提出了DPD模拟方法，其运动方程中力的部分类似于布朗动力学，包含耗散力和随机力。Hoogerbrugge等人认为分子动力学在粒子运动上过于细致和繁琐，因此DPD方法中的粒子不再对应真实的粒子，而是对应流体的微小区域或高分子的几个链段。在此基础上，Espanol等人进一步得出随机力和耗散力权函郑州大学硕士学位论文数必须满足耗散涨落定理，从而形成了目前的DPD方法。 JonesJL等人[’9]利用耗散粒子动力学模拟了液滴受到剪切力作用时附着在固体表面时的变形情况，并适当地模拟了固液界面的粘附边界条件，以及在一定的剪切流条件下液滴脱离固体边界的情况。同时，还考察了液滴在固体表面的接触角等问题。 clarkAT等人[20]利用耗散粒子动力学对重力作用下悬垂的单个液滴以及受剪切力作用时的液滴变形及破碎的情况进行了耗散粒子动力学模拟。重力作用下悬垂液滴的变形情况与通过求解LaPlace方程所得结果能够很好地符合;受到剪切力作用时的液滴变形情况与几ylor理论相符。Torza等人[2‘]利用DPD模型进行试验，结果与理论上液滴破裂的整个过程十分相似。陈硕等人1221论述了耗散粒子动力学在复杂流体流动方面的研究进展，并简要的介绍了并行计算技术在耗散粒子动力学方面的应用。 冯剑等人123】在Groot和W价en耗散粒子动力学的改进veriet速度算法基础上，通过优化方法获得了一个简单的计算兄的表达式，进而提出了一种自洽获取该值的优化改进Verlet速度算法。陈硕、赵钧等人1241利用耗散粒子动力学对三维微通道中液滴的动力学进行了研究。通过提高两种流体之间的守恒力参数，使这两种流体互不相溶，并分析了这两种流体在Flory一Huggins理论中相应的参数。此外，目前针对微通道内生物巨分子悬浮物的耗散粒子动力学研究取得了初步的成功，Fan、Phan一Thien和ehen等人[2，J引入有限拉伸非线性弹性(F即E)的珠簧链模型和蠕虫模型(Worm一like)对生物DNA分子悬浮物在微通道内部的流动进行了研究，并得到DNA分子在突扩突缩狭缝通道中多种典型的构形。由于耗散粒子动力学本质上属于离散粒子动力学方法，它要求巨大的计算量才能展现真实流体的物理特性，一个有效的解决办法就是采用高性能并行计算技术。针对耗散粒子动力学所描述的系统，可以采用并行技术来实现对较大系统的真实模拟。和被分配到各个微处理器上的粒子数量相比，参与通信的粒子数量是很小的一部分，因此可以保持较高的并行效率。所以开展耗散粒子动力学的高性能并行计算研究，也是目前耗散粒子动力学研究的一个方向。总之，耗散粒子动力学作为一门介观尺度的计算技术，模拟研究工作仍然处于起步阶段。郑州大学硕士学位论文5课题主要内容首先，对DPD方法的基础理论进行系统的学习，掌握了DPD方法的使用原理，并且吸纳了一些关于DPD方法改进的观点;其次，对简单流体进行基于DPD方法的模拟验证DPD方法在流体模拟领域应用的可行性。再次，对聚合物流体进行基于DPD方法的模拟，得出关于聚合物流体性质的一些结论。最后，对聚合物流体的流变性能进行分析，证明了采用DPD模拟的方法也可以很好的测定聚合物流体的流变性能，结合经典流体力学方程和DPD理论方法介绍了聚合物流体的模拟过程的介观和宏观藕合分析方法。.6本文结构第一章绪论部分介绍不同数值模拟方法，分析DPD方法目前研究情况，明确课题的研究意义，确定了主要研究内容，形成一个比较清晰的结构脉络。在第二章中，我们会详细介绍这种新兴的数值模拟方法—耗散粒子动力学方法，分析其在整个计算模拟领域的地位和研究范围，给出它的理论依据，陈述相关的一些算法以及数据处理方法。对简单流体的流动模拟会在第三章中叙述，通过调整合理的参数，与己有结果进行对比，验证程序的有效性，并进行讨论。第四章进行聚合物流体基于DPD方法的模拟，介绍如何利用DPD方法模拟测定聚合物流体的流变性能，将模拟结果代入流体力学经典方程中进行宏观一介观祸合分析得出一定结论。最后本课题做出总结和展望，耗散粒子动力学用于聚合物流体的研究尚处于开始阶段，新兴的耗散粒子动力学的应用会有很大前景。郑州大学硕士学位论文第.二章 耗散粒子动力学理论基础耗散粒子动力学(DissipativepartieleDynamies，DpD)是以统计力学为基础的介观尺度流体模拟方法。在DPD系统中，基本单元是被称作粒子(P叭icle)的离散动量载体，它们在连续的空间内运动，并采用离散的时间步长。这些粒子是一些具有质量的粗粒化实体，每颗粒子的运动代表大量分子的集体动力行为，这样研究就可以在粗粒化粒子尺度上进行，而忽略更小尺度的结构和运动情况。DPD方法的主旨在于解决宏观尺度的传统方法运用到复杂流体系统时所出现的困难，同时，又改善了微观尺度模拟方法，如分子动力学方法，所需要的计算量巨大的弱点。本章先列出耗散粒子动力学方法的基本模型，后介绍各种常用的处理方法，最后对模拟过程进行概括总述，进行程序结构的描述。 2.1DPD基本理论2.1.1标准DPD方法DPD方法的基本模型，采用N颗粒子充满整个连续空间的DPD系统来展开叙述。每一颗DPD粒子由坐标气，动量p，和质量m，来描述。DPD粒子的运动由牛顿运动方程式描述为，对于任意一颗DPD粒子i，表达为:de dy:=v;，一“dt__。(2.1)一dt=f;+厂“其中，弓和v，为粒子i的位置和速度矢量;厂为所有其他粒子(不包括粒子本身)作用在粒子i上的力;F“为粒子受到的由聚合物链模型弹簧杆提供的作用力，t为时间。DPD系统中假定所有的粒子都有同样的质量，采用无量纲单位质量作为粒子的质量，即m，=m=l。郑州大学硕士学位论文1镶图2.1粒子之间的相互作用范围 FigZ.lTheinteraetionrangebetweenPartieles在DPD系统中，每颗粒子在一定范围内和周围粒子发生相互作用(如图2.1所示)，它们之间的相互作用力(即公式2.1中的厂)由保守力、耗散力和随机力三部分组成，表达式为: 厂·艺岭十十丫)(2.2) 其中歼为保守力(eonservativeforee)·衅表示耗散力(dissipativeforee)、叮表示随机力(randomforee)，其中角标扩表示粒子i与粒子j之间的两两相互作用。式中的求和符号表示对一定的截断半径rc内所有的不包括粒子i的其他粒子进行求和(粒子之间的距离几>r.，则粒子之间的相互作用为零)o截断半径被定义为无量纲单位长度，即r.=保守力耳是一种软排斥力，作用在两粒子质心连线方向上，表达式为:.‘了、/、|l之1、了.||anU刀︸︺ F 1一鱼r，份
:时，wj)一0，咋为粒子之间相对速度，气=v，一vj，:为粒子之间相对速度为耗散系数，:前耗散力的作用总是阻止粒子之间的相对运动，因此耗散力作用的结果会降低系统的动能。系统动能的降低将由随机力叮引起的粒子的随机运动进行补偿，随机力可表达为: 丫=，“认从r.;力权重函数，当‘>几时，w“=小易为随机变量，有以下性质:(2.6)式中:二为随机力系数;护认)也是与粒子之间距离。有关的权重函数，称为随机(易卿·0，(易(t肠〔’》=伍*戈，十她、歼一才’)且i护j，k护l，随机力大小是时间相关的积分 [e认wa犷=只丫()·t叮(+:牺t(2.8)与保守力一样，随机力和耗散力也作用在粒子质心的连线方向。在前人的研究中提到，粒子之间的动力相互作用是由耗散力和随机力来描述，这两部分作用相互补偿从而使得系统的平均动能为一恒定值。为了满足扩散波动定理，使系统保持平衡，要求式郑州大学硕士学位论文 (2.4)、(2.5)中的两个权重函数之间、耗散力系数厂和随机力系数，之间相互关联，关联式中的变量只能任选，一个，称为系统的平衡条件。关系式表达如下: 研日=[n认w犷犷=2成。T夕」〕(2.9) (2.10) 则有 式(2.10)中:K，T为系统的Boltzmann温度，将K。T作为能量的无量纲单位K。T二犷二2厂(2.11)权重函数并没有指定的表达式，但是无论何时都必须为正，比较典型的形式为:，份表示全体取平均。等式中第一项表示由于DPD粒子的动能传递形成的力，第二项表示粒子间相互作用力。压力就由压力张量的迹得到:1 _(2.16)二一一才r万 2.1.2改进DPD方法Fan等人124]提出了一种改进的DPD方法。主要的目的就是提高DPD模拟流体的粘性及斯坦顿数。但是这种办法会大幅度增大DPD方法的计算量，其方法如下:在DPD理论中，耗散力与随机力权重函数要满足下式:(l一、丫，、<: 式中，rc之1，:‘2，当:=2时，为标准DPD二次权重函数。由耗散理论我们知道，DPD系统的粘性受两部分的影响，一是粒子的扩散系数，二w”认’一卜“认‘·{oPD.十，今之rc(2.17)是粒子间的耗散力。受耗散力影响的粘性粉”起主要作用。DPD流体的粘性和斯坦顿数(Schmidt~ ber，Sc)的计算公式如下:11一—+—PD1一—1+—+粉D二——1、}(2.18)s+15+25+35+4s+5/ (2.19)式中，p为DPD流体密度，D为扩散系数。郑州大学硕士学位论文Fan等人[24l对DPD粘性公式分析得出::越小，r.越大，则粉”越大。当:=2时，，’低丫一瓜。。当‘=05时.，”低丫一2%6’树5增灯78倍.。同时，增加rc也可以增加粘度。但是增加rc会大幅增加每个粒子与周围粒子之间的作用力。而所需要花的计算时间正比于截断半径的立方，因此这样做会增加模拟的计算时间。虽然增加耗散力强度z也可以增加粘度，但是通过实际的模拟检验，这种做法会增加稳定模 拟所需温度的困难[24】。Fan等人提出将，减少到0.5，凡增加到1.5，会极大的缩短DPD模拟的反应时间，并且会使流体的粘度及斯坦顿数增大数倍。2.2模拟基本原理在陈述耗散粒子动力学方法的模拟过程之前，先介绍一些模拟中常用基本设定。这样处理可以有效优化算法，使我们的模拟更加简便，模拟效率提高。2.2.1无量纲方法在计算机模拟中通常都采用无量纲单位来表现我们所研究的问题，进行这样处理的优点如下:A.数据处理时方便，减少误差。由于模拟涉及到的都是极小的量级，如微米级，直接应用实际值工作时十分繁琐，工作量巨大，还会造成很大误差。采用无量纲单位，就可以保证数据的一致性，数据简单，处理方便;B.简化方程。方程中的一些参数，在经过无量纲处理后就成为单位参数，在方程中可以将其隐去，如此处理可使方程得到简化;C.防止由于计算机硬件问题产生的计算溢出。对于微尺度问题的数据，都是0.0000001到0.000000001之间的数级来表示，对于一般的处理器很容易出现溢出，这样会使结果受到很大影响。使用无量纲方法可以有效的避免溢出现象。耗散粒子动力学模拟中所有的物理量我们均同时进行无量纲处理。根据参考文献[z6一271中的叙述，对所研究的问题进行无量纲单位与实际单位的转换(以水为例)，从而将模拟问题与实际情况联系起来。用X%表示参数的实际单位量，比如t%表示实际时郑州大学硕士学位论文间:用[x1表示参数J￥的单位量，比如卜1表示单位长度;〔1表示单位时间等;用X表示无量纲量，如p表示无量纲密度。三者之间的关系为:x%=x[x}(2.20)例如截断半径的实际长度为:几%=1.608、10一6。，在DPD模拟中通常令截断半径为单位长度，即几=1，由rc%一rc卜1=1.608、10一m，推导出单位长度为:[z]=1.608、10一6m，根据单位长度可以方便地将无量纲长度与实际长度任意转换。很多参数的换算关系可以由已知基本量纲间接计算得到。按照上例同样的方法可以推导出所有参数的无量纲单位与实际单位之间的关系，我们将其整理成表格，便于比较。在耗散粒子动力学的无量纲单位转换中，基本单位是长度卜}、密度巨1和能量卜]，其他单位都是由基本单位推算得出;然而传统的基本无量纲单位是长度卜!、质量卜1和时间卜l。需要强调的是耗散粒子动力学模拟时，各个粒子并不代表真正的原子或分子，它代表流体体系中一小部分区域。因此在耗散粒子动力学模拟时，DPD单位采用的是无量纲化后的单位，可以通过选择标定单位将真实的物理量对应到DPD体系中。此时真实的位移、速度和时间(;，，，t)可以由DPD模拟体系中的(;，，，t)来表示:一r一V一t尹一冈，’一伺，‘一厕(2.21)式中，[rJ，[vj，[tj为选用的位移、速度、时间标定单位。在模拟结束后，可用模拟的数值乘以标定单位得到对应的真实数值。DPD模拟结果就是用粒子的质量、位移、速度、温度、相互作用参数等物理量表示的。 郑州大学硕士学位论文 表2.1DPD方法中各参数无量纲单位与实际单位的换算关系基基本参数数实际单位无量纲单位单位量 长长度rc%=1.608x10一mmrc=ll{月=rco/0crc===1.608x10一16mmm 密密度(水))户%=1.oxlo，kg/m，P=44!p]=p%/ppp===0.25x10，kg/m， 能能量:%=凡T=4.115、10一，‘egrg((T(=298K)))￡==K。T==ll刀刀!:}·内/““=4.115x10一14e馆馆 质质量【m}=【户』{Vlll===l·040x10一”kggg 速速度时时间[v]动面向向==3.45x10一3m/ss 【，1=【Ll!v』』==4.66lxlo一455根据表中所列参数的无量纲单位换算关系，我们把研究问题无量纲化;相反的，如果要知道一个无量纲模型所对应的实际问题情况，只要知道无量纲单位换算关系，就可以得到实际问题的尺度。2.2.2聚合物FENE型结构高分子是含有原子数很多、相对分子质量很高、分子很长的巨型分子。高分子由很大数目(10，一l护数量级)的结构单元组成，每一结构单元相当于一个小分子，这些结构单元可以是一种(均聚物)，也可以是几种(共聚物)，它们以共价键相连接，形成各种结构。研究高分子结构的意义在于了解分子内部和分子间相互作用的本质，即通过对于分子运动的理解建立结构与性能间的内在联系，掌握结构与性能的关系。单就分子中所含原子个数、相对分子质量的大小和分子的长短还不足以表达高分子的结构特性。人们通过长期的实践和研究，证明高分子是链状结构。每个高分子里含有郑州大学硕卜学位论文面面函一1.面一-.面一面面..函面..‘面..亩面面面面‘‘.面‘‘百面面‘面面面面面一面一面~面声‘面丽‘‘可声‘面亩‘一‘~‘面面一个或数种原子或原子团，这些原子或原子团按照一定的方式排列，首先是排列成很多重复结构的小单元，称之为结构单元，再通过化学链连成一个高分子。如果高分子的分子链呈不规则的线状(或者团状)，高分子是一根根分子链组成的，就称为线型聚合物;如图所示。这类聚合物分子间没有化学键结合，分子长链蜷曲成团，或者伸展成直线，这取决于分子本身的柔顺性及外部条件。在受热或受力情况下分子间相互移动(流动)，因此这类聚合物可以在适当溶剂中溶解，加热时可以熔融，易于加工成型，且这种特性在聚合物成型前后都存在，因此可以反复成型，习惯上也称这种材料为热塑性塑料。图2.2(a)直链线型(b)带有支链线型(c)体型 Fig2.2(a)linearPolymer(b)branehedPolymer(c)bodyPolymer 在DPD方法中，我们将聚合物用珠一簧链结构(FENE结构模型)来表示，如图:翻目确O闷丫州日以R 图2.3FENE珠一簧链结构示意图Fig2.3F阴Ebead一springehaineonfiguration珠簧链包含M个珠子和戈二M一1个具有一定长度的弹簧.珠子代表熔体中相互作用的质点，没有重量的弹簧代表链中的约束.其中，弓(i二1，2，3，…，M)表示珠子的位形向量，g(，.=I，2，3，…，从)表示珠一珠之间的连接向量，r..表示珠簧链的质量中心，尺、表示珠簧郑州大学硕士学位论文链的轨迹半径。DPD方法中，FENE模型中的珠子由DPD方法中的粒子(particle)代替，则作用在粒子i与粒子j之间的作用力可为: 可=-门1.氏1.2又.火、/、.1.(2.22)其中，H是弹簧常数，几是粒子i与粒子j间长度，‘是最长链段长度。我们在时间上定义两个时间常数，一个是弹性珠簧链的时间常数，另外一个是刚性珠簧链的时间常数，分别为:杏 代12K，T(2.23)ha=~万丁丁，勺=峥厂1这两个时间常数的比值一般我们用来表示弹簧的有限拉伸常数: U :_3几_K叫兄HBT二二—二二—(2.24)如果我们用共来表示聚合物链段的平衡长度，则聚合物总的链段平衡长度为恤。一，妈，对于FENE弹簧链结构来说， L。=，区2.2.3粒子初始位置Vb+5(2.25)耗散粒子动力学方法对流体的模拟采用的是将粒子按照接近真实物质结构的方式排列，由此得到初始状态的所有粒子位置。常采用面心立方体(FCC)格式和体心立方体(BCC)格式。面心立方体格式就是指粒子在模拟区域内按照面心立方体的格式排列，即在每个立方体的顶点和六个面的中心分别各有一颗粒子(如图2，4a所示)。体心立方体格式就是指粒子按照体心立方体的结构排列，每个顶点和体心上都分别有一颗粒子(见图2.3b)。在我们的模拟区域内采用面心立方体(FCC)的结构排列，所有粒子的初始位置便可以确定。郑州大学硕士学位论文 {{、声矛， 卜卜、矿.. {{，才少夕Vt、、\\一、目目‘三_一___玉_一~__图2.4(a)面心立方体结构示意图; 伪)体心立方体结构示意图Fig2.4(a)Theconfigu时ionof伽eeentercubie;助Theeonfigurationofbodyeentereubie2.2.4粒子初始速度对于流体粒子，一般来讲速度的给定是随机的，但是要满足以下三个条件:第一个条件是，动量必须守恒，即:1~丝{ 衅“俨一方菩俨第二个条件是，动能与热力学温度有关，即满足下式:(2.26) 合客二酬式中，a代表x，y，z三个方向。刀兀BT2(2.27)第三个条件是一个选择条件，即速度初值要符合Max叭lell一。】忱mann分布。 P沙少二!二一二丁二二1exPI一二二亡言I‘、f。，、奋「l。，v?二，1戈乙从Bl/L艺人Bl，」(2.28)M~n一Boltsmann分布给出了质量为的原子i在温度T下沿x方向速度为vi，:的概率。为了使粒子的初始速度满足该分布，可以根据(式227)先求出平均速度。则每个粒子的初始速度为平均速度乘以该粒子的Max叭吧11一。h刀刀越m分布概率。郑州大学硕士学位论文2.2.5边界条件及处理方法在模拟中，由于计算机能力的，所能处理的体系是有限的。也就是说，模拟必然是在一个有限的区域内实施，而实际区域相对于计算区域来说是无限大的。模拟中粒子会碰撞边界并尝试着逃离模拟区域，区域边界附近的粒子就会增多，而实际系统只有一小部分粒子在壁面附近，这就是所谓的“表面效应”。为消除模拟区域 引入的表面效应，需要引入与分子动力学类似的周期性边界条件(PeriedicBoundaryConditions，PBC)，这样就可以构造出一个准无限大的体积来更代表宏观系统。 占占OOO00占OOO006OOO00奋... 黔黔O叹乡叹乡睁III翰0厂、01.给O..二))) 石石OOO00吞OOOOO石OOO00奋... 称称00，招O奋...弓O 、奋奋香香OOO!言”奋OO移O一下奋勺勺香‘公OOOOOOOOO0书 0000公0书称0公图2.5二维周期性边界条件体系 Fig2.5ThesystemofthetwodimensionPeriodbound哪eonditions以二维周期性边界条件为例:如图2.5中，粗线的格子为模拟格子，应用了周期性边界条件后，消除了界面效应，但模拟结果和体系真实的性质之间仍然存在误差。当模拟区域内部存在较大体积时，模拟结果会受到模拟区域的影响，这就是“有限尺寸效应”为解决这一问题，理想做法是让这个体系可以随着时间的变化根据需要实时改变模拟区域的大小，但这样使模拟变得复杂，而且给分析结果带来困难，因此这种做法不可行，要采用固定的边界条件。一个选定的体系是否会受到尺寸效应的影响而无法得到可靠的结果，这需要进行合理的验证，具体方法是:保持模拟区域的体积不变，改变其长，宽，高的比例，通过观21郑州大学硕一}:学位论文察这个体系的压强张量的变化来衡量体系中的有限尺寸效应，如:对于一个完全处于平衡态的体系，它的压强张量的对角线上的三个元素应该是相等的，即:sxx二sy;=又z[2l8。如果S二、Syy、522不相等，就说明选定的模拟区域存在尺寸效应，就要对所选择区域进行调整。常见的模拟中除了这种只要在三个维度方向实施标准的周期性边界条件就可以实现的无限大系统以外，还有一种就是存在壁面的情况。这种系统的模拟只能在两个维度方向上实施周期性边界条件，壁面需要采用特殊的方法实现。最初采用“冻住”的粒子(“冻住”指在OK下，所有的粒子没有振动动能和位势能)实现壁面的固定不动，这些粒子可以与自由粒子相互作用，但是位置固定在原处没有动能。这样的方法有三个问题:其一，由于耗散粒子动力学方法中的粒子采用的是“软”作用力，没有分子动力学方法中的约束力那么强，自由流体粒子就会穿透壁面，造成流体粒子的流失;其二，壁面附近运动的粒子与壁面之间会产生速度的滑移，同时密度会产生波动;其三，在靠近壁面处温度会降低，因为冻结的壁面会使附近的流体粒子冷却。由于这些原因，关于壁面边界的处理一直是一个难题，先后出现了很多种解决方案。Kon梦9]采取增加壁面粒子密度的方法解决粒子的穿透问题，但这样会使排斥力增大，导致壁面附近粒子减少，模拟不能反映出真实情况，且密度波动严重。之后Jones提出了改进方法，可以保持壁面和流体密度相同，在壁面附近又建立一层强排斥力层，从而增大壁面与自由粒子排斥力的方法，这样做虽然可以防止粒子的穿透，但会导致接近壁面处的密度扭曲，而且也没能实现无滑移边界条件。Revengalsol提出可以产生反射的固体壁面粒子层，当相互作用力不足以阻止粒子穿透壁面时，反射层就会发挥作用将粒子反射回流体内部。他的论文中又提到三种反射方式:镜面反射，反弹反射和麦克斯韦反射。在镜面反射的前后，壁面切线方向的速度分量不变，仅仅是法线方向的分量变为反向;在反弹反射中，切向分量、法向分量都变为反向;而麦克斯韦反射是指反射回到流体内的粒子速度要服从麦克斯韦分布。Revenga分析了影响壁面滑移速度的主要参数，在DPD系统中五个关键的无量纲参数是:m(粒子质量)、:(耗散力系数或摩擦系数)、rc(截断半径)、凡T(温度)和丫=。一肠(粒子间的平均距离，其中d是空间维数，n是密度)。无量纲摩擦系数的定义如下:郑州大学硕士学位论文 :二二dvT(2.29)很小的情况，第一种和第三种反射方式会产生很明显的速度滑移，而反弹式却能表现良好。既便如此波动。V:=厚指热速度。当:很大时，三种反射方式都能得到无滑移边界，但对于:它们的温度曲线在接近壁面处仍然都会出现异常，同时密度产生强烈的willemsenl3’}采取在模拟区域外额外增加一层粒子的方法，延续壁面边界，建立合理的速度分布。组成这一粒子层的粒子的位置和速度由临近壁面截断半径内的一层流体粒子来确定。例如，为了实现固体壁面零速度，模拟区域的附加层粒子的切向、法向速度与壁面内部的一层粒子反向。当DPD粒子撞到固体边界，就会产生反弹反射，这种方法在忽略保守力的时候表现很好，但是当考虑保守力时，就会出现密度波动。第二层(附加层)设置在边界外rc和Zrc之间，利用这一层产生排斥力，防止粒子穿透的同时改善壁面附近的密度波动。相对来说，这种采用两层边界层实现固体避免的方法有一定的效率，但还是存在着些许边界滑移。Duong一Hon梦21提出的方法与willemsen的方法相似，在模拟区域外增加两层粒子，如图2.6，当粒子穿过壁面，就会产生新的速度嵘，二2呱I，一琉J(气aI，表示壁面的速度矢量)反弹回模拟区域。粒子层与边界之间的距离为a、rc和azrc(0)，斌BTl一V、了、了、了、+一，j(2.32)PK，:·争’工了(()krr2dr郑州人学硕】:学位论文为了得到状态方程，Groot令密度p为1一8之间以0.5的间隔递增的一系列数，排斥力系数选择三组数据:a=15、a二25、a=30，进行计算。由此得到的数据经过处理和分析后发现，当p>2时，三种系统成为同一条曲线，经过数据拟和可以得到压力的近似公式: 尸·PKoT+an，a(a=0.101士0.0001)将此公式代入公式2.31，得到K一，=1一Zaap/K。T、1+0.2ap/K。T(2.33)对于己知的水的压缩系数K一，二16，推导出ap/KBT二75。原则上密度的选择是自由的，但是与每颗粒子相互作用的粒子数随密度的增加线性增加，每一步所需要的CPU方增加而递增，由于效率的原因，一般选择尽可能小的密度值。时间随密度的平在Groo‘的论文中选择。·，，一2，、;对于其他密度一75凡%。根据这个原则，我们在模拟中采用。二18.75的数据。根据参考文献[32]壁面的排斥力系数。*=5.0，壁面粒子与流体粒子的排斥力系数a，=扣苏可=”w.““。算法中涉及到的参数和模型参数，基本上都按照上述方法选定。2.3.3系统的运行确定了模拟参数之后就可进行正式模拟计算。先在不施加外力的情况下让系统自由松弛，逐渐达到平衡状态，达到平衡态所需要的时间称为弛豫时间。当系统达到平衡状态以后，进一步施加外部干扰，进行非平衡态计算。运算中最重要的就是计算2.1中讲述的模型方程，用数值求解的方法解微分方程组，采用Verlet方法，方程组2.1变成:。(t+动=2。()一t。(t一叫十伽)，at)(其中，a(t)即辉仓)/m，由于质量m采用无量纲数，设定为1，这样方程表达为:布王件月tJ，‘f一△l这两步的空间坐标位置及t时刻的作用力，就可以算出下一步t+Al时刻的坐标位置。为了将其写成更加简洁的形式，我们令 ‘，·n△‘，。(”)·rr(t，)，月(·)=只(n)t(2.35)郑州大学硕士学位论文这样方程组简化为: :(二，)=2:(·)一。(一，)+月(·)(A，)，，(i=l，2，…，N)由空间坐标又可以算出粒子的运动速度为(2.36)已知一组初始空间位置御0，工分’好，则通过求解方程组一步步得到认(，’王认(3’}，…。 v)·，=认‘二”一。‘一，’)/2△才这里在第n+l步算出的速度是前一时刻，即第n步的速度。因而动能的计算比势能的计算落后一步。(2.37)由以上方程组设计得到的模拟步骤如下:(1)给定初始空间位置价‘“，{认‘，’好，(i=1，2，…，，)(2)计算在第n步时粒子所受的力脚·，}(3)利用公式:沙·‘)=2n()一nn(n一，)+砂)伽丫，计算在第。+1步时所有粒子所处的空间位置妙·”}(4，计算第n步的速度:讨·，·仓‘二，’一。‘一，’)/2△，(5)返回到步骤(2)，开始下一步的模拟计算由以上分析可知，要真正求出微分方程组的解，需要给出两组空间位置俨，}和认。’}，而且动能和势能的计算不同步，这样，最好采用改进的Velocity一veriet方法，可以修正上述两个问题。公式如下: :(二’)一:(·)+△tv;(·)+生△才)，:(·)2(’Vi(2.38)一(n+l) =讨”)+;△城(·)一(n+1)i(2.39) 月(·‘，)=汽[:(”‘’)，V(2.40) r”J-1、(.、l__‘_、_‘_.:、_v户”’‘产=v全“，+一△t}户;“‘夕+户;“‘，‘，}2根据此算法得到的模拟步骤为:(2.41)郑州大学硕士学位论文(，)给定初始空间位置认‘。，}，(i=1，2，…，、)(:)给定初始速度俨}(，，利用公“:。(”十‘，=:(n，·△。(n，+告协t)’;，计算“第n+‘步”所有粒子所”的空间位置认‘二，’}(“)计算在第n+，步所有粒子的速度何n+l 讨一’=讨·，+合八‘。:‘·，·:‘”‘”，(5)返回到步骤(3)，开始第n+2步的模拟计算。度，并且数值计算的稳定性也提高了。按照以上步骤，设定足够多的时间步，就可以得到一组模拟结果。(2.42)这种修正的算法比前一种的算法好些。不仅可以得到同一时间步上的空间位置和速2.3.4宏观物理量的计算实际计算宏观物理量是沿着相空间轨迹求平均来计算得到的。例如对于一个宏观物理量“，它的测量值应当为平均值万。对于初始位置和动量为奋‘“’(0}和夕“)，(0}(上)标N表示系综N个粒子的对应坐标和动量参数)，求解具有初值问题的运动方程，便得到相空间轨迹({r‘“’以御“’(t)})。对轨迹平均的宏观物理量，的表示为 万=吻尚恻{r‘“’以份‘“，介))}如果宏观物理量为动能，它的平均为(2.43) 瓦·吻尚卜‘·护，哪在时间的各个间断点产上计算动能的平均值(2.44)由于在模拟过程中计算出的动能值是在不连续的路径上的值，因此公式应该表示为 、=六豁嘿以此为例，宏观物理量可以采用积分求极限的方法计算求出。(2.45)郑州人学硕士学位论文模拟最终目的是研究流体的流动特性，因此得到的数据都是与流体特性相关的各种参数，如速度、压力、温度、密度等。经过通过与已有理论和实验解进行对比，调整算法和程序，不断完善模拟方法。2.3.5模拟程序程序开发:编程思路包括三部分，分别为:初始部分、平衡部分和数据采集部分。A.初始部分的任务a.定义所有变量.，给定模拟所需的必要参数，如各个力学公式中的参数、熔体和管壁的密度、模拟温度、各种模拟步长、模拟区域、边界条件等。并根据设定的参数计算出其它相应的参数，如模拟所需的粒子数。b.给定所有粒子(包括流体、管壁、聚合物)的初始位移。c.给定所有粒子(包括流体、管壁、聚合物)的初始速度。d.使用提高效率的方法。e.计算最初的力。B.平衡部分的任务a.计算粒子新的速度和位移。b.计算每个粒子受到的力。c.矫正粒子速度。d.运用边界条件。“.标定速度。f.达到某个任务步长时，执行这个任务。9.输出体系状态，如程序已执行的模拟步数、模拟体系的温度、动量等。h.检测是否达到平衡步数，达到则进入数据采集部分。没有达到平衡步数则返回平衡部分任务a，继续下一个循环。郑州人学硕」学位论文C.数据采集部分的任务a.达到采集数据的步长时，保存数据值。b.计算粒子新的速度和位移。c.计算每个粒子受到的力。d.矫正粒子速度。e.运用边界条件。f达到某个任务步长时，执行这个任务。9.输出体系状态。h.是否达到模拟总步数，如果达到模拟总步数，处理保存的数据，然后输出数据处理结果，最后结束程序。没有达到则返回数据采集部分任务a，继续下一个循环。郑州大学硕士学位论文程序运行流程图:程序大致分为两部分:将粒子随即分布于模拟区域，得到所有粒子的初始位置;计算粒子之间的相互作用和运动速度及产生位移，实现对实际问题的模拟。流程图如下:开始输输入模拟区域大小小 壁壁面粒子排列 建建立主流区主流区粒子排列生成粒子位置文件输出结果图2.8模拟区域程序流程图 F192.8SimulationofregionalProeess郑州大学硕十学位论文开始读读入粒子位置文件件输输入初始数据设定模拟参数条条件判断平平衡判断每每一时间步计算生生成各参数结果文文件输输出结果图2.9土程序流程图F192.9ThemainProgramProeess 郑州大学硕士学位论文2.4本章小结通过DPD方法基本原理的介绍，各种处理方法的分析以及模拟过程的叙述，勾勒出了耗散粒子动力学方法模拟的整个轮廓。DPD粒子反映一个小区域内的位置和动量，而模型没有给出粒子与粒子间的明确边界以及粒子内部分子原子之间的情况，而是忽略这些细节当做粗粒来处理。对于这种模型，为满足动量守恒的条件提出耗散力衅，即描述粒子边界的相邻区域之间由于相对速度产生的粘性力。对于实际流体，衅会导致流体区域内的升温，而对于DPD流体，这种能量可以看做是能量从粒子外部这种宏观模式流入粒子内部的微观模式，这种效果可看作能量被移出系统。随机力丫可以理解为由于内部温度的原因导致的流体随机波动。在DPD流体中，表示能量从微观模式转化为宏观模式。保守力叮表示流体区域之间的压力表现，对于真实流体，这种压力由微观尺度上粒子之间的相互作用和随机运动产生。在DPD流体中，这种力被认为是粒子之间的一种排斥力，可以分析出:耗散力描述了能量从宏观模式流入微观模式，在效果上使系统降温;随机力描述能量从微观模式流入宏观模式，使系统温度升高。根据DPD平衡理论，系统必须保持平衡状态，那么对于两个具有相同的总动量而每个粒子动量不同的系统，平衡时两系统在任意方向上的运动概率相同就要求耗散力和随机力相互关联达到这种平衡，即叽认)=叫(ur)的关系式。郑州大学硕士学位论文第三章简单流体流动模拟耗散粒子动力学方法对于流体流动方面的研究还很不够，缺乏DPD方法与已有的经典分析解以及传统数值模拟解进行比较，因此我们应用DPD方法对流体的动力学行为进行了研究。3.1泊肃叶(Poiseuille)流模拟我们采取的模拟模型是流体在两无限大平行平板间的层流流动，应用DPD方法来实现这种基本的流动形式，并与己有结果[34】进行比较，验证程序的正确性。3.1.1泊肃叶(Poiseuille)流/夕牛//////////////// 图3、1poiseuille流示意图Fig3.lThePoiseuilleflow平板间的粘性层流动运动虽然简单，却能够揭示粘性流动的一些本质。如图所示，两平板管壁之间充满了粘性系数为声的流体，在流场中施加一恒定的压力梯度，则流场速度的解析解为:v=一丽丽(万一了少式中，，为速度，p为压力。h，咖‘::，、‘3·1，郑州大学硕_卜学位论文3.1.2模拟方案泊肃叶流动模型是在Z方向上壁面存在，X、Y方向是无限大的粘性流动，采用前文所提出的三层粒子和无滑移边界条件处理Z方向的壁面，在X、Y方向实施周期性边界条件，采用简单的DPD粒子模拟粘性牛顿流体的流动。基本方案确定后，首要问题就是确定合适的模拟区域。在相关论文[3，’241中采用的模拟区域的大小为:60x3x30，这样选择的结果是符合要求的，可以较好地描述研究的模型。由于我们所模拟的结果要与己有结论进行对比，检验程序的正确性，因此也采用与论文中相同的情况: 一30兰x<30，一1.5兰y<1.5，一15三z<巧。在模拟区域内将粒子按照FCC的格式进行排布，由此可以确定所有粒子的初始位置。粒子的总数取决于模拟区域的大小、几何形状以及流体和壁面的物质密度:流体的初始速度是根据给定的温度随机设置，而壁面采用“冻结”的粒子，初速度为零。整个系统共采用11800颗简单DPD粒子，其中10000颗充满整个流动区域，1800颗为壁面粒子，壁面分别由平行于xy平面的三层粒子组成，最内层位置为:=士15。接下来选择合适的参数，每一个参数的确定都按照2.3.2中介绍的方法，模拟中涉及到的所有参数见表3.1:表3.1重要参数目录 模模拟区域壁壁面粒子数数60X3X3000n=180000 流流体粒子数数时时间步长外外力单单位质量截截断半径单单位能量n。=1000000△t=0.0222g=0.0222阴=111rc=lllKBT11=1 密密度P二333 随随机力系数数a二3郑州大学硕_l学位论文 VVeriVet算法系数流流体排斥力系数壁壁面排斥力系数兄=0.6555a/=18.7555“w二5.000 流流体与壁面间排斥力系数数平平衡时步长数a可=9.6888n。。=200000 计计算.总步长数n;。a、=2.0x10555系统的运行使用Velocity一Veriet算法，计算每一步的位置、速度和加速度;计算粒子间相互作用力时采用2.2.6中描述的元胞分割和关联表法。在统计数据时，将流体区域的Z方向分成300个小格，将数据按照格子为单元进行记录，描述流体特性的数据是将每个单元内每隔时间步长10‘的样本数据进行统计平均，对这些数据进行处理后就可分析出系统的流动特性。郑州大学硕一1学位论文3数据处理我们将计算得到的数据文件进行处理，得到了流体性质的各种表现。A.速度分布曲线 40——t二1000&2000t二400t二800————t二2000t二4000t二10000——一t二80003.5 ——t二20000t二200 >x2.0a二18.75，仃二3， g==0.02，p二3 一15一10一5051015图3.2泊肃叶流动在Z方向的速度发展曲线 F193.2Thede、‘eloPmentofveloeityProfi!esinPoiseuilleflowofasimPleDPDfluid模拟得到的速度发展曲线如图3.2所示，处理数据时，我们每隔200个时间单位观察一次，最终选择了最能表现发展过程的曲线组合在一起。根据图中速度的发展趋势可知，在1000时间单位后流动到达充分发展阶段。对于充分发展段，我们选择t=2000时的曲线与理论解进行对比。牛顿流体的充分发展速度分布理式为:「 11}2/_、2]1}v:=vmax}l一f一11[\\n夕」(3.2)共丁极禺阴一丰，_v_max=_v_，WJ乙J_甲.h_即。_带_于，_.__距。__*。、.，品少=一二一一~户动，=工。二。二。匕声日甲足候戮_卜__趴_，、_特l、_拟上.、，、lL1_据与理论结果进行对比，发现两者相互吻合。郑州人学硕卜学位论文U勺一曰4口勺 :n只︺—--3.525*【1一(2115)**215imU!ationdata︼/ 2.0飞考\\\\a二18.75，a二3 1.0粼仔冷z0.5g二0.02，p二3、飞﹄ 夕一15一10双议认叭 0.0图3.3充分发展速度分布曲线与N一S方程的理论解比较 F193.3Then」llydeveloPedvelocityProfileeomPareswelltotheN一5PredietionsB.压力分布a二18.75，a二3 40一g二0.02，p二3dg巴jSn一15一10图3.4泊肃叶流动的压力图 F193一落ThePressureProfilesofasimPleDPDfluidinPoiseuilleflow郑州人学硕_卜学位论文 压力是根据2.1中的公式(2.15、16)计算得到，在完全发展阶段的压力分布如图3.4所示，同样在主流区表现良好，在靠近壁面处有一定的波动情况，这种不稳定不会影响到主流区。C.剪切应力分布剪切应力的分析解是一户岁，将DPD流体模拟所得的计算解与理论解的进行比较(图3.5)，可以看到除了接近壁面处的一些偏离，其他主流区吻合得较好。 ——5imU!stionn 1.SJ一pgz 0.5类(D。刀一0.5 一1.5 一2.015一10一5051015Z图3，5泊肃叫流动的剪切应力分布 F193.5ShearstressdistributionofasimPleDPDfluidinPoiseuiileflow我们的模拟得到了描述泊肃叶流动的各项参数表现，也与理论解进行了比较，发现结果完全吻合。郑州大学硕十学位论文3.1.4结果讨论 最后我们将得到的结果与Duong一Hong在论文中的计算结果进行比较，如图3.6所示，其中3.6(a)是Duang一Hong的模拟结果，3.6(b)是我们的结果，3.6(c)是Fan的结果。2二Ve!oe扮P动Ie.V:Slmol.t朗一v厅以(1·(乍5)2)弓510Zo1015图3.6(a)Duong一Hong模拟的充分发展速度分布图 F193.6(a)ThefollydeveloPedveloeityandtheNavier一Stokessolution ————3.525*!1一(Zj15)**21115imU!ationdataaaZ图3.6(b)本文模拟的充分发展速度分布图 I了193.6(b)ThefullydeveloPedveloeityinoursimulation4l郑州人学硕卜学位论文曲一翻漏~漏一一森一~~一磊一山~动‘叶愁.了5.叮二3力、 万多J麟诱，一‘内丈月丫违呼曰，愧1‘月︸4透奋硕礴孟确书莽t圣L.滚峨!1，l禅‘l宁.吞240O2t鑫26(犯钾粗‘卜粉r叨﹄f爹咨卜t—翻m川每娜d翻怕·“拍甲钾低25九.冬O石O一立琦仍Z图3.6(c)Fan模拟的充分发展速度分布图 F193.6(e)ThedeveloPmentofveloeityProfi]esofasimPleDPDfluid比较图3.6(a)与图3.6(b)，二者的模拟参数设定相同，施加的外力大小不同:Duong一Hong的模拟g=0.01;本文选择g=0.02。因此可知速度发展趋势是相同的，只 有速度大小不同。在Duong一Hong的论文中，最大速度vmax二1.756;本文模拟结果v，a、n一3.5280根据Duong一Hong的模拟结果，当参数g与我们的参数也相同时，推算出最大速度v，nax=3.516，两者误差很小。说明我本文的模拟结果是正确的，从而验证了程序的可行性。与Fan论文中的结果进行对比，二者的模拟区域大小相同，所有的基本参数选择相同，施加外力也相同，不同之处在于Fan的权重函数w。(r)=w吴(r)一(l一:)，，而我们选择了w)·r(，压，“前文初步交代了权函数的不同选择是可以满足不同的要求，后者比前者产生更大的耗散力。在这里通过对比分析不同的权函数可以实现怎样的不同需求。Fan的速度最大值vma、=8.639;本文的模拟中vmax=l，759。由此可知，平方根的权函数可以模拟比平方的权函数粘度大的流体，也就是说，模拟粘度大的流体时，如果平方的权函数不能满足要求，就可以采用平方根的权函数。从理论角度分析，平方权函数小于平方根权函数，前者计算得到的耗散力比后者小，这就是前面提出的平方根权函数可以得到更大的耗散力的原因。耗散力是一种描述粒子郑州大学硕士学位论文边界的相邻区域之间由于相对速度产生的粘性力，也就是说耗散力大即粒子之间的粘性力大，因此，采用平方根权函数模拟得到的流体粘度大。通过与理论分析解，也充分验证了程序的可行性。3.2方形通道内流动的模拟利用DPD方法模拟方形通道内部流体的充分发展流动，在介观尺度上再现连续介质流体力学的研究，然后用基于N一S方程的数值模拟结果验证模拟的正确性。3.2.1模拟方案方形通道尺寸为:一巧‘x<巧，一10兰y<10，一10兰z<10。通道在X方向为流动方向，计算中在该方向实施周期性边界条件，流体在Y和Z方向受壁面约束。每个壁面依然由三层粒子组成，在壁面附近实施了无滑移边界条件[34l。模型总共采用了60000个DPD粒子，其中48000个为流体粒子，其余的12000个为壁面粒子。计算区域在x、Y和Z方向共设置lx5Ox5O个格子，时间步长设置为0.01，重力加速度g二0.01在X方向施加给每一颗流体粒子来驱动流体流动。其他参数的选择与3.1.2中表3.1相同。系统运行时同样采用Velocity一Veriet算法计算每一步的位置、速度情况，用元胞分割和关联表法计算粒子间相互作用力。3.2.2结果讨论选取充分发展段的速度分布曲线与应用Flueni软件模拟得到的情况进行对比，图7方形通道内充分发展流动速度分布图比较。郑州大学硕十学位论文「IUentDPD图3.7方形通道中充分发展速度分布 F193.7ThefullydeveloPedve!ocityProfileinasquareehannel由图可见，DPD模拟结果与Fluent软件相比所得到的数值模拟结果符合很好，这说明DPD在模拟由简单粒子所表示的流体方面与经典的流体理论相符合。DPD所模拟的方形通道内充分发展流体流动结果表明，DPD用于简单流体的计算中，能够较好地再现经典流体力学所描述的流体行为。这也表明了DPD方法可以在介观尺度范围对流动问题进行模拟研究。3.3本章小结本章模拟了简单流体泊肃口十流动和方形通道内流体的三维充分发展流动，并与理论分析解、经典文献结果进行对比，验证了算法的正确性;对方形通道流动的模拟证明了DPD方法可以实现对连续介质流体力学的模拟研究。郑州大学硕一卜学位论文第四章聚合物流体流动模拟聚合物流体的流动行为十分复杂，并且直接影响到高分子材料加工工艺的选择及高分子材料使用性能的发挥。因此，在聚合物成型加工工作中需要研究聚合物流体在外场作用下的流动行为。描述聚合物流体流动行为的表征数据有:粘度，剪切应力，切变速率等等，而这几种数据均与聚合物流体的剪切流动关系密切，此外，我们在聚合物加工过程中常常出现因剪切作用出现残余应力，从而影响制品质量。因此，我们本章就对聚合物流体的剪切流动进行模拟，并就一些相关问题进行进一步的讨论。4.1剪切流动模拟为了实现剪切流动，我们采用的模型是两无限长平行平板反向运动形成定常的剪切流动，采用这种流动形式能够产生稳定的剪切流。 4.1.1模拟方案模型平行平板采用“冻结粒子”构成，聚合物流体用简单的DPD粒子来实现，在X、Y方向实施周期性边界条件，在Z方向为了防止粒子穿透壁面以及速度滑移，}司样采用三层粒子边界，令两平行平板以一定的速度向相反方向移动形成稳定剪切层流。对于模拟区域，我们选择一20兰x<20，一10‘夕<10，一10‘:<10。在模拟区域内将粒子按照FCC的格式进行排布，由此可以确定所有粒子的初始位置。粒子的总数取决于模拟区域的大小、几何形状以及流体和壁面的物质密度:流体的初始速度是根据给定的温度随机设置，而壁面粒子初速度为零。整个系统共采用11800颗粒子，其中10000颗充满整个流动区域，1800颗为壁面粒子。模拟中所采用的参数与表3.1相同，流体密度为p=3;随机力系数a=3;耗散力75;壁面与流体间排斥系数厂二4.5;壁面排斥力系数a，二5.0;流体排斥力系数价二18.力系数a叮一9.68;时间步长△I=0.02;Veloeity一Verlet系数兄=0.65。对于剪切流动，自然要提到描述剪切流特性的一个重要参数—剪切率。它的大小郑州人学硕士学位论文与平板速度有关，定义式为: D=卫三兰三h(4.1)其中，v:表示上板速度，v一:表示下板速度，h表示两板间距离。设定了平板的速度，根据这个公式就可以得到剪切率。在本中，上下板速度大小相同，方向相反，速度暂设定为0.5，即v，二一v_，二0.5。参数都设定后进入正式的模拟。系统的运行依然使用Velocity一Veriet算法，计算每一步的位置、速度;计算粒子间相互作用力时采用元胞分割和关联表法。4.1.2模拟结果及数据处理模拟得到的速度分布曲线如图4.1所示，处理数据时，我们每隔200个时间单位观察一次，发现速度发展很快，t=200时就已经稳定。由于采用了Fan所提出的无滑移边界条件，在靠近壁面边界处没有明显的速度滑移，也不影响主流区的流动。流动性质的其他参数曲线见图4.2、4.3，其中图4.2中为剪切应力的分布情况;图4.3是压力图。 ——t二2000 ——t二4000t=t10000>X住 1.5十一---，--~--一了一---一一10一50510Z图4.1剪切流动速度分布图 F194.1ThevelocityProfilesinCouerteflow郑州人学硕卜学位论文 S利的S﹂)Qa二18.75，a二3，Vx二0.5，p二3SXZ{一1十一一 ·10一一--了------尸----一丁一一一一一一50-.一下------尸一一-一~〕510图4.2剪切流动中剪切应力分布曲线 F194.2TheshearstressProfilesinCouetteflow日曰”Il曰日以nll，户︿dg二SJjgS︸︹︿、卜日钊“丫任二18.75，a二3， Vx二0.5，p二3O十一 一10rses一e一-leses一5，w一e-一，一0---1一5一~-州10图4.3剪切流动中压力分布曲线 F194.3ThePressureProfilesinCouerteflow郑州大学硕士学位论文经过观察并结合模拟结果进行分析，我们确定所模拟的剪切流动是稳定成功的，为我们接下来的进一步研究讨论提供良好的基础。 4.2FENE链构象变化在剪切流动中，聚合物粒子组成的FENE链将随时间逐步拉伸并伴随有自由旋转。为了观察FENE链在流动中，随时间所发生的拉伸、自由旋转和取向等变化，我们在计算机上利用Matlab、’叭sualC++软件对它进行了可视化处理。图4.4、4.5分别为链在不同时刻位形变化图，图4.4为单根分子链，图4.5为100根链。户(a)t=0才， (b)t=0.02郑州大学硕卜学位论文才.争 (e)t=0.04图4.4 一根分子链不同时刻位形图F19.4.4Theconfigurationof1moleculeehain(a)t=O气沪萝(矛.20=t)b厂郑州大学硕士学位论文了声尸令声 (e)t=0.04 图4.5100根分子链位形变化图F19.4.5Theeonfigurationsof100moleeuleehains图4.4一4.5所示为链的构象随时间变化，链在流动方向上的移动在三维坐标下显示，从图中可看出分子链在移动过程中发生了拉伸、自由旋转等变化。4.3聚合物流变性能聚合物的流变性是指在外力场作用下，溶液粘度与流速或压差之间的关系。就象固体聚合物的力学行为可以用模量(定义为应力和应变之比)描述一样，聚合物的流动行为可以用粘度(定义为应力与应变速率之比)表征。普通溶液的粘度是材料的一个常数，它只与温度和压力有关，而与形变速率及时间无关。但聚合物溶液的情况复杂得多，粘度随形变条件而变化，而且聚合物溶液的流动伴随着弹性效应，因为施加于体系的部分能量以可回复能量的形式被储藏起来，因此粘度还与时间及形变速率有关，即聚合物溶液是粘弹性的。4.3.1聚合物的粘性通过聚合物流体的剪切模拟，我们得到了比较理想的模拟结论，通过模拟的结果我们进一步进行分析，希望在此基础上对于一些高聚物流变性能进行讨论，如粘度等。为了得到高聚物流体的粘度，一般情况下我们大都是通过采用实验手段，通过使用流体进行仪器对待测定流体进行实验测定得到结果，比如通过使用旋转粘度仪，落球粘郑州大学硕士学位论文度计等。但是实验测定粘度的过程常出现因一些外界因素的干扰从而导致粘度测定结果的不准确，因此我们一直希望找到更为理想的粘度测定方法，尽量避免出现测定误差。本文在此方面也进行尝试，即:通过流动模拟的方法从微观角度来测定粘度，以期得到更为准确的粘度。在本文中，前面介绍过应力张量可以根据Irving一Kirkwood模型进行计算: 百一万(今阴Z“‘·“‘刀+1/书N﹂了，t.P﹃，卜艺j心1>、Jl、、zla心几刀(4.2)其中，V是系统体积，。，是DPD粒子质量(m，一1)，N;是粒子数目，从。和与粒 子，的特定速度，例如:u，。=v心一可(x)，石(x)是x处的流体流动速度(streamveloei‘y)，<二>表示全体取平均。等式中第一项表示由于DPD粒子的动能传递形成的力，第二项表示粒子间相互作用力。从上式中我们得到应力张量，因此，我们使用本构方程可以将粘度表示为:S二 粉=D(4.3)式中，凡:表示剪切应力张量，D表示剪切率，在上一节当中已做介绍。我们通过模拟得到Sxz，通过模拟参数的设定进行计算可以得到D，这样我们就可以利用上式准确的计算出流体粘度，从而避开了利用实验手段测定所可能带来的外界影响和结果误差。4.3.2聚合物的弹‘陛高聚物在剪切流场中，除有粘性外，还表现出奇异的弹性，存在法向应力差。根据第一、第二法向应力差私、NZ定义:第一法向应力:N，二S二一凡:第二法向应力:NZ一又:一Syy我们根据模拟可以通过计算得到S二，S。，S二，因此聚合物流体的法向应力差也可以通过模拟得到。郑州大学硕士学位论文由此我们可以推出:第一法向应力系数:第二法向应力系数:少l=普梦2=告其中树l和梦2分别定义为第一法向应力系数和第二法向应力系数。通常笋2二一0.1梦，，因此，粘度函数加上第一、第二法向应力差系数就可以完整描述稳态剪切流场中的粘弹性流体的应力状态。4.4宏观一介观尺度的祸合4.4.1流体力学基本方程A.连续性方程连续性方程(continutiyequatino)是质量守恒原理在流体运动中的表现形式。张量形式可表示为: 零、v·伽)·”口t(4.4)式中p为密度;y为速度矢量。B.动量方程 动量方程(momentumequation)是动量守恒原理在流体运动中的表现形式。剪切流动的粘性流体动量方程:Dv，___ P一万午=j，一VP+V·几口t(4.5)右边第一项表示体积力，第二项为压力，第三项为应力。用此动量方程与连续性方程对于一定的流体类型，在具体的边界和初始条件下，可以求解一些高聚物熔体加工中的流动行为，并可求得速度、压力、应变和应变速率等流变参数的变化关系。郑州大学硕士学位论文4.4.2宏观一微观拱合分析上一节中我们介绍了流体力学里面两个基本方程:连续性方程和动量方程。我们在前面的工作中利用DPD方法对聚合物流体的流动过程进行了模拟，并且可以根据DPD方法取得一定的模拟结果，作为主要的就是可以知道粒子的速度，位移，以及流体的应力等，根据所得这些结果我们可以计算流体的粘度。对于连续性方程和动量方程我们采用张量方式进行表达，如下: 助.r7了__一、n连续性万程:二厂十v‘又尸v少=v口t(4.6)Dv， 动量方程:P不万Ij丁一几一助+v’公，(4.7)在DPD方法中，应力张量表示为:l 万=一一艺心i艺、l、召z、矛la几介刀(4.8)通过模拟我们可以求得剪切应力张量Sx:，因此，压力就由应力张量的迹得到:l _(4.9)P=一一tr百将所得剪切应力张量S二和压力p代入动量方程中，则连续方程与动量方程组成二元方程组，这样我们就可以对速度和总的外力两个未知数进行求解。对于聚合物FENE模型，其弹簧力F“在前面我们已经介绍过为: 可一亡在DPD方法中，外力的和为:(4.10) f=关+F，(4.11)其中，厂为所有其他粒子(不包括粒子本身)作用在粒子i上的力;厂为粒子受到的由聚合物链模型弹簧杆提供的作用力。DPD粒子的运动由牛顿运动方程式描述，对于任意一颗DPD粒子i，表达为:郑州大学硕士学位论文 dv:__。一-二=f;十户’Jdt式中，科和从为粒子i的位置和速度矢量。(4.12)通过我们对连续方程与动量方程组成二元方程组进行求解的结果，我们可以得到粒子的速度的大小，因此我们可以由速度求出剪切速率从而得到分子链作用力，再由分子链作用力求出新的应力，将其带回DPD方法中，这样我们又可以得到速度的大小，知道了速度的大小我们又可以利用DPD模拟得到一个新的应力值，这样一个循环过程就实现了介观一宏观两个尺度上的祸合计算。图4.6宏观一微观藕合分析路线 Fig.4.6CouPledMacro一MieroSimulation郑州大学硕】学位论文4.4.3棍合结果讨论上一节我们介绍了宏观一微观祸合的具体思路和计算方法，现在我们具体将祸合方法进行应用，对进行和没有进行祸合分析所得的模拟结果进行对比，如下图: ，峥n月自U门三月1l.门nAveragestressvstime0”UJ八、.一一日!一一““n.，n.卜.旧.LI1 nocouPledeouPled乏︶叻祠戈5﹄︵减﹄。只自乙f︸，月Un0八︸︶”0 0l.we. 0.0090.014爪me/s0，024图4.7祸合效果对比 F19.4.7CouPledeffeeteontrast通过对比我们不难发现，采用祸合分析的结果比没有采用祸合的结果更能准确的反应聚合物流体流动过程中的的应力变化。因此，我们采用祸合分析的方法能够很好的体现出流体流动的真实情况，让我们掌握更准确的信息。4.5本章小结在本章中对聚合物流体进行剪切流动模拟，得到了流体性质的结果。通过整个过程的实践，对使用耗散粒子动力学方法模拟聚合物流体的流动过程有了更加深入的认识。我们通过模拟结果对于聚合物的流变性能进行了一定程度的分析讨论，证明了采用DPD模拟的方法也可以很好的测定聚合物流体的流变性能，还可以通过模拟结果得到聚合物FENE链的构象的具体变化。郑州大学硕士学位论文我们结合经典流体力学方程和DPD理论方法介绍了聚合物流体的模拟过程的介观和宏观祸合分析方法，说明DPD方法在流体问题的研究上具有很大的应用前景。郑州大学硕士学位论文第五章总结与展望本文的主要工作是应用耗散粒子动力学这种介观尺度模拟方法，对聚合物流体的流动进行模拟，大致分成三部分来实现:首先，介绍应用耗散粒子动力学方法进行模拟的过程。其中详细介绍了DPD基本模型，模拟中常常用到的处理方法，如何确定初始条件，选择合适的参数，以及最终数据的统计平均方法，描述研究问题的物理量计算等。其次，应用本文方法模拟了简单流体泊肃叶流动和方形通道内流体的三维充分发展流动，并与理论分析解、经典文献结果、Fluent软件的模拟结果进行对比，验证了算法的币确性。再次，实现对聚合物流体的剪切流动进行模拟。并且通过模拟结果对于聚合物的流变性能进行了一定程度的分析讨论，证明了采用DPD模拟的方法也可以较好的测定聚合物流体的流变性能。最后，我们结合经典流体力学方程和DPD理论方法介绍了聚合物流体的模拟过程的介观和宏观祸合分析方法，说明DPD方法在流体问题的研究上具有很大的应用前景。经过一系列的研究，实现了DPD方法在聚合物流体流动的模拟，同时体现出它在处理流体模拟的独特优势。但是，对于不同组分的聚合物流体以及多组分的聚合物流体的流动如何运用DPD方法进行数值模拟我们还需进一步研究‘能得到准确的解决方法。此外，如果较好的解决了上述问题，那么我们可以根据研究结果进行基于DPD方法的聚合物流体数值模拟软件开发工作，这也是有一定意义的一个可深入研究的的工作方向。总之，对于采用DPD方法进行流体流动模拟研究尚处于初级阶段，还有很大的空{。习需要去探索。郑州人学硕士学位论文参考文献 lrr..乙勺-FloryP.J，JAmChemSoc，1936，58:1977.YOrk，1953.FloryP.J，PrineiPle，;ofPolymerChemistry，CornellUnivPress，New 【3」.StoekmayerW.H，J一ChemPhys，1952，9:69. !4」.FloryP.J，JAmChemsoe，1941，63:3083.!51，StoekmayerW.H，JChemPhys，1943，11:45.!6」.MaeoskoC.W，MillerD.R，Maeromolecules，1976，9:119 17」.NakanishiH.，PhysRev，1980，B22:2466.!8」.FellerW，AnIntroduetiontoProbabilityTheoryanditsAPPlieations，JohnWiley&Sons，NewYork，1957. 【9」.M.P.Allen，D.J.l’ildesley，ComPuterSimulationofLiquids，1987{1Ol.N.MetroPolis，A.W.Rosenbluth，M.N.Rosenbluth，ete.SegmentationsofSPatio一TemPoralImagesbySPatio一TemPoralMarkovRandomFieldModel，JChemPhysl953，21:1087. l川.W.EvanGunsteren，H.J.C.Berendsen，AngewChemIniEdEngl，1990，29:992!12」.D.C.RaPaPort，ThepreSS，Cambridge，1995.ArtofMoleeularDynamiesSimulation，CambridgeUniversity 【131.GCieeottiandW.G.Hoover，Molecular-DynamicsSimulationofStatistiealMeehanicsSystems，North一Holland，Amsterdam，1986 {141.D.FrenkelandB.Smit，UnderstandingMoleeularSimulation:FromAlgorithlnstoAPPlieations，Aeade汀一iePress，London，1996. 115」.U.Friseh，B.HasslaeherandY.Pomeau，Phys.Rev.Lett.1986，56:1505. t161.R，Benzi，S.SueeiandM.Vergassola，Phy.ReP.1992，222:145.【17」.S.Sueci，R.BenziandF.Massaioli，Int.J.Mod.Phys.C，1993，4:409.{18〕.P.川oogerbruggeandJ.M.VA.Koelman，SimulatingmieroseoPiehydrodynarniephenomenawithdiSS;ipativepartieledynamies，Europhys.Lett.，1992，19:155一160. 【19」.JonesJL，Lal入1，RuddoekJNandSPenleyNA，DynamiesofadroPataliquid/solidinterfaeeinsimpleShearfields:amesoseopiesimulationstud又1999，112:129一242. 【ZOJ.ClarkAT，LalM，RuddoekJNandWarrenPB，MesoseoPicsimulationofdroPsin郑州大学硕士学位论文致谢本论文是在我的导师申长雨、曹伟两位教授的亲切关怀和悉心指导下完成的，他们严肃的科学态度，严谨的治学精神，精益求精的工作作风，深深地感染和激励着我。从课题的选择到论文的最终完成，二位教授都始终给予我细心的指导和不懈的支持。_三年来，他们不仅在学业上给我以精心指导，同时还在思想、生活上给我以无微不至的关怀，在此谨向二位老师致以诚挚的谢意和崇高的敬意。感谢国家橡塑模具中心其他老师给予我的支持和帮助。感谢国家橡塑模具中心810室的张羽标、岳鹏、李刚、胡德富、张亚涛同学和室友陶沈健、沈永龙、马国发同学以及模具中心06级的其他同学，三年的朝夕相处中，我们相互学习、彼此关怀，建立起了深厚的友谊，度过了人生中最美好的一段时光，愿我们的友谊地久天长!最后还要感谢培养我长大含辛茹苦的父母和一直支持我的朋友们，谢谢你们!张文羽中二零零九年五月 