CPU-GPU 耦合的多尺度模拟应用.doc

CPU-GPU 耦合的多尺度模拟应用摘 要:文章简单明了的阐述了多尺度问题以及研究手段。论述了常见多尺度问题的模拟计算方法与研究进展，对现有研究的局限性和存在的问题进行分析，指出了进一步研究多尺度模拟与计算的必要性。离散颗粒系统（DPM）已经被普遍的用于颗粒流体系统的相关模拟当中，然而DPM在超大型系统展开模拟的时候出现计算量特别大的问题，这一点极大的阻碍了它的全面推广应用。不过这些年来初露锋芒的图形处理器（GPU），由于其具备极强的计算能力以及不错的并行性，给离散颗粒系统的推广应用开辟了另外一跳道路。关键词：多尺度；模拟；DPM；应用1第1章 引 言1 颗粒材料的概念1.1定义与分类颗粒材料，英文名称为granular materials, granular matter，学术界的定义是通过众多相互接触在一起的颗粒所构成的系统，因为在低于1微米的情况下，热运动极其明显，因此颗粒尺度往往会现在大于1微米的范围里。在超过1微米的情况下，热运动即可不用考虑在内，颗粒与颗粒之间产生的相互作用力仅仅只有接触作用力，换言之，量子力学在这种情况下可以忽视，仅仅只要考虑经典力学领域的计算。按照有没有存在液体一般可以把颗粒材料区别成干颗粒材料以及湿颗粒材料两种类型；按照结构组分又可以把颗粒材料区别成单相以及多相颗粒材料两种类型；这里面的提到的单相颗粒材料，具体来说就是拥有单一屋里性质的相关颗粒构成的干颗粒类型的集合体，而多相颗粒则是多种物理性质不一的干颗粒构成的集合体（或者称之为多相颗粒流），以及具备间隙的流体的颗粒组合体。在非理想状态下，因为颗粒与颗粒之间往往都有空气充斥其中，那么从这一点来说，单相颗粒材料只存在于理论当中，在现实生活中是不可能有的。不过要使研究更容易开展，通常会人为的对间隙当中的空气不作考虑，默认当前颗粒材料属于单相颗粒材料，这种情况一般出现在空气影响比较轻微的时候，包括干砂粒堆以及小麦堆等。颗粒材料根据力学行为的不同又可以区别成静态（static)和动态(dynamic)，而动态一般又包括三种小类：其一为准静态(quasi-static)，其二为动态(dynamic)，其三为快速流(rapid flow)。按照固体颗粒浓度的不同以及孔隙的差异还可以把颗粒材料区别成密相颗粒材料(dense granular materials )以及松散颗粒材料(loose granular materials)还有稀薄颗粒材料(dilute granular materials )三种类型。另外，依照颗粒浓度以及无量纲的剪切速率还可以把颗粒流区别成准静态、过渡态与快速态这三种类型。1.2 颗粒中力的传递处于相互接触、尺寸大小一致、形状相同的均匀小颗粒共同构成的颗粒组合体内部形成的有序的、对称的二维层面的颗粒组合体，所谓加压实验，就是把颗粒组合体放进具备水平地面的刚性结构当中，利用从垂直方向上给予集中压力的方式，详细观察并记录颗粒组合体内部相关力的扩散情况。在随机理论的基础上，构建力的传递模型。这套模型在结构上来说是二维有序以及对称分布的。在这套模型当中，所有颗粒都受到了从上方相互接触的两个颗粒所传递的力，同时这个力又通过利用相同的途径往下再次到达与之接触的邻近颗粒，具体可以看图1，处于处的颗粒1把受到的力往下传递到达邻近的点上的颗粒2以及处的颗粒3。图1颗粒1受力并传递到达下层颗粒2与3按照能量传递的理论，可以得出相应的传导一扩散方程： （1）这里面，分别表示波动量基于方向上的偏移系数以及扩散系数。假设代表力产生的首个部分，那么转移概率密度可以表示为 （2）方程2代表的含义是作用在最开始处即的力波动量最终可以传递到处颗粒上的概率大小。力的波动量基于方向上的平均位移可以表示成，同时存在有的Gaussian涨落。由此不难发现，当力在颗粒堆积体系不断往下传递的过程中，力波动量基于方向上的传递会由一个常数来决定，传递到平均位置的偏差系数同样伴随颗粒深度的不断提升而增大。假设，即对称分布，偏移系数，此时扩散方程就变成 （3）因为颗粒组合体当中存在搭拱效应，同时力波动量基于方向产生偏移转移，导致集中施加的垂直压力在传递过程中伴随着颗粒深度不断提升而不断向两边扩散，在一定深度的颗粒堆积体内部，会发现在集中力作用线的位置，也就是说颗粒组合体的中心部分，不是受力最大位置。综上所述，在进行定性分析的时候，这套理论模型能够与实验相匹配。要注意一点的是，由块状颗粒构成的颗粒组合体在内部力的传递上不符合连续介质弹性理论的基本条件，也就是说，其相对于连续介质在力学特性上存在区别。因为离散特性的存在，粗粒料的多数相关现象都无法通过用已有的连续体力学模型来进行解释说明，因此连续介质模型中的变量在离散颗粒系统当中的任何位置都不具备连续的特性，只具有统计意义。1.2 颗粒材料相关模型1.2.1 颗粒材料宏观层面的连续体模型所谓连续体模型普遍的适用范围在密相颗粒材料，一般包括经典连续体模型以及Cosserat(微极，micro-polar)连续体模型这两大类型。颗粒材料的连续体描述已经在岩土工程和化工领域中得到广泛采用。 (a)经典连续体模型。经典连续体模型已经普遍运用到化工行业以及土木工程领域。最早在1650年的时候就对土壤颗粒展开了相关研究，法国人Coulomb在这个层面发表了土的抗剪强度以及土压力滑动模型，这种简单的模型之后加入到散体极限破坏相关的莫尔一库仑(Mohr-Coulomb)标准。在这以后，在经典连续体模型不断发展的前提下，土力学逐渐发展并形成了一门新的学科领域，其通常用于土壤颗粒，还有间隙流体相关的研究工作。一直到现在为止，土力学已经发展出了一系列数学模型以及相应的本构模型，自然土力学相对应的理论以及数据模拟研究也不断发展并前进。不过，美中不足的是，经典连续体模型并未设定内部长度参数，也就是说无法体现微观结构层面的信息，也就不能对剪切带进行预测。一旦出现应变局部化，同时在应对相关问题，包括静态建筑维护方面、慢速变形层面中出现的土壤蠕变以及地基问题，还有快速流层面的滑坡问题的时候就会面临尴尬的局面。在化工行业当中，涉及的颗粒材料往往具备流体性质，而以类Navier-Stokes方程最为基础的颗粒材料流体系统普遍使用到类似流化床等一系列领域内。 (b)微结构连续体模型。从微观层面来看，各种结构属性，包括颗粒排列情况、相互接触情况，这些都跟颗粒材料在宏观层面的力学理论存在密不可分的关系。考虑微观结构影响有两种途径，即离散途径(discrete approach)与微结构连续体途径(microstructural continuum approach)。微结构连续体途径主要是利用宏观角度的连续变量来转化成宏观角度上的材料本构属性，而这里的宏观连续变量一定要可以体现出微观离散性质。相当多的学者通过这一途径来探寻颗粒材料的基本特性，比较典型的当属Cosserat连续体模型。该模型增加了旋转自由度的概念，同时在此基础上得到与之对应的偶应力，还有相应的特征长度，而通常来说，特征长度就可以通过本构方程一定程度的体现一部分微观结构信息，具体来说如颗粒尺寸信息等，因此如果利用该模型来模拟颗粒材料内部的微结构是最好不过了。不过经典Cosserat连续体模型同样有着不足之处，这里面一个明显的问题就是这套理论模型仅仅只可以用来部分的体现微结构信息，对于排列信息等则显得无能为力，同时其和颗粒材料离散结构特征之间并不能形成很好的显示关联。 Nemat-Nasse等研究者加入了结构张量(Fabric tessor)来修正之前的连续体理论模型，该模型的关键是如何建立颗粒材料离散的微观力学与宏观上连续行为之间的联系。利用结构张量这一属性来体现颗粒材料微结构特性，这样就可以很好的解决之间存在的局限性，并且在宏观应力和结构张量之间搭建了一座桥梁，帮助本构关系可以很好的体现颗粒材料微观层面上的结构变化过程，也就是在这个前提下，颗粒材料基于微观力学连续途径研究领域取得了长足的发展。这条途径利用应力以及应变两个概念诠释了颗粒材料在宏观层面发生的变形行为，在颗粒与颗粒之间存在的局部运动及其相关力学行为的基础上，形成了微观本构理论关系。20实际90年代，Chang针对颗粒材料进行了划分，具体分割成表征元以及微单元与接触三种层面，这三个层面分别对应宏观应力应变情况、局部应力应变情况、力和位移的基础关系。假设将颗粒材料基于微观层面的离散特性加进来，那么之前提到的经典连续体相关的应力应变定义就不成立了，正是在这个背景下，更多的研究集中在怎样描述颗粒材料的微观结构信息上，还有怎样全面的定义微观结构上的应力应变与选择哪一种均一的形式来搭建微观颗粒特征基于宏观行为两者间的联系。Chang CS等在Vogit均一化理论的基础上应用最优拟合假设(best fit hypothesis )、运动假设(kinematic hypothesis )、静态假设(static hypothesis)等一系列理论模型给颗粒材料在宏观以及微观两个层面上的联系，还有相应本构关系展开了一系列基础研究。Combou M具体的阐述了用于颗粒材料均一化方面的各种经典途径。Kruyt NP, Satake M, Bagi K等人则推出了颗粒材料相关各种几何说明以及应力应变体系。 微结构连续体理论体系属于广义上的连续体系统。Rothenburg等人经过研究表示，如果将颗粒的旋转情况排除在外的话，那么就可以搭建仅仅只有相关颗粒尺寸信息的一套宏观上的连续体体系，而如果要考虑到颗粒旋转这个因素，那么宏观层面上就会搭建微极类(micro-polar)以及Cosserat类这两种连续体系。而Chang CS为代表的相关学派则利用对颗粒微观层面上的位移信息进行处理，包括颗粒产生的平动以及转动等，进一步组建宏观位移场，在这个基础上得到一种能够把离散系统等效的变更为连续体系统的途径。在这个过程中，等效的连续体系统可以很好的体现离散系统的全部特征，同时，还能够退化变成微观层面下的连续体系统的全部相关类型，如高阶梯度微极连续体(high-gradient micro-polar continua，考虑颗粒的自旋及旋转)、高阶梯度偶应力理论一Cosserat连续体(high-gradient couple stress，加入颗粒旋转情况)，还有高阶梯度非微极连续体(high-gradient non-polar continua，排除颗粒平动与旋转情况)等。Jansen在19世纪末期针对沙漏以及粮仓底部平面压力进行的研究就率先发表了颗粒材料在连续途径上的描述办法，其人为的定义颗粒材无论在哪种平面的情况下，其内应力都处于均匀分布的状态，同时认为垂直向产生的应力相对于水平向产生的应力存在正比的关系。而到20世纪末，Bouchaud等学者旨在描述沙堆表面产生的崩塌现象针对性的创建了BCRE系统，其核心内容就是把系统当中的沙堆表面看成是两个类型，其一是稳定的颗粒层，其二是滚动的颗粒层，这两者又能够在位错(dislocation)以及粘合(sticking )的形式完成互相转换的过程。静态颗粒要产生位错，一般要依靠重力来进行，也可以是运动颗粒之间的相互碰撞产生，在这个同时，一旦颗粒运动到相应的安息角以下的层面，那么颗粒又会变成稳定状态。在这以后，更多的研究人员在此系统的帮助下，进一步的诠释了颗粒材料表层流动现象、颗粒分离与分层现象的基本原理。而为了将颗粒离散属性纳入模型当中，Copper Smith在1996年率先创造了后来被称作q-model唯象的概念。1.2.2 颗粒材料的微观离散系统Cundall以及Struck利用牛顿经典力学理论，围绕接触力展开计算，创造了离散颗粒体系，也就是熟知的离散单元办法，这套系统已经广泛应用于颗粒材料相关行为的微观尺度研究领域。最开始颗粒材料离散体系包含了下列几大假设: ( 1 )颗粒为圆形且视为刚体;(2)颗粒与颗粒之间可以在接触点层面产生一定量的重叠；(3)颗粒与颗粒之间产生接触面积认为是无限小的;(4)颗粒与颗粒之间的接触点在重叠量大小的确定上认为和接触力线性存在关联，前提是颗粒没有发生分离或者活动行为；(5)颗粒间滑动遵循Mohr-Coulomb定律。这套理论系统又可以称作离散单元形式的软球系统。在这之前Campbell创建了硬球系统，同时将该系统应用到剪切流的分析研究上。Kishino在进行颗粒材料准静态的模拟过程中，创建了颗粒单元方案(granular element method)。而Thornton则在接触力学的基础上针对离散单元接触本构及其程序展开了进一步的发展与扩充。离散单元系统型的关键就在于接触模型。应用比较广泛的就是弹簧-粘壶-滑片模型，是通过Cundall等人创建。由于离散单元系统的不断丰富与发展，加上大范围的投入使用，越来越多人将研究方向定位于接触模型及其具体算法上。这里面比较显著的成果就是Oda以及Bardet发现的颗粒材料在剪切强度还有剪切带的发生发展都和颗粒发生滚动行为高度相关，而Iwashita以及Oda在Cundall创建的接触系统的框架下又添加了经过修正以后的离散单元系统，同时在该系统中添加了滚动摩擦阻矩这一概念，同时又考虑了颗粒之间的滚动机制。Vu-Quo等人则对接触力的发生发展远离进行了详尽的研究，尤其提到一点，要注意接触面的相应半径相对于接触力大小之间量化关系的确定，在这个基础上又得到了基于弹塑性变相关的接触本构关系。而Feng等人则重点提到了滚动摩擦力基于颗粒物理行为相关模拟研究中的关键程度，并且进一步说明颗粒和颗粒之间发生摩擦行为的物理原理直到现在在学术领域依旧缺乏得到统一认可的系统模型。Jiang等人则对纯滚动以及纯滑动行为展开了重新定义分析研究，其认为基于接触点发生的位移能够分解得到相应滚动位移以及滑动位移，通过数学模型详细的推到得到Iwashita等为代表的修正接触系统，还有将接触宽度产生的影响涵盖在内的重新调整之后的离散单元接触本构系统。另一方面，因为Cundall是在光滑接触模型的基础上创建得到离散系统，那么Fortin等为代表的学者就以非光滑接触模型为背景创建了对应的非光滑接触离散单元系统，希望通过此举来增强高接触力在计算方面的精确度。虽然说在学术界已经广泛认为离散单元法可以帮助科研人员在微观层面进行颗粒材料行为的相关行为研究，把颗粒材料结构物的整个宏观域模型化为单个颗粒的集合体进行数值模拟，哪怕是在当今超级计算机飞速发展的背景下，也是无法实现的。那么科学的应对办法就是将微观层面上的离散颗粒系统和宏观层面上的连续体系统通过某种中间形式有机的结合到一起，在这种背景下，颗粒材料的多尺度模拟系统就诞生了。1.2.3 分子动力学模型(Molecular dynamics,MD) 20世纪50年代末期，Alder以及Wainwright率先通过分子动力学硬球系统理论对气体以及液体相关状态方程展开研究，在世界上首次从分子动力学的层面来研究物质宏观特性。在这以后，分子动力学理论在物理、化学还有材料科学等学科上发展上都产生了深远的影响。分子动力学基于颗粒材料模拟中的应用普遍集中在颗粒速度比较高的高能颗粒系统当中。根据处理接触方式上的区别，分子动力学系统通常包括两种风格各异的体系，其一是在时间步驱动下开展的TD体系(Time-step driven MD)，其二是在事件驱动下的ED体系(Event driven MD )。从根本上而言，分子动力学系统在颗粒材料领域的应用相对于颗粒离散单元系统是一致的，尤其是这里面的TD体系在实际操作环节上和离散单元可以说是完全一致。而ED体系则是在建立在颗粒与颗粒之间仅仅只有瞬间碰撞的基础上，利用相关的碰撞系统来得到颗粒的动量变化情况。由于计算机系统及其技术的不断进步，分子动力学系统与并行技术有机结合到一起已经成为颗粒材料模拟领域的下一个热点。2 多尺度模拟概念作为化学工程的基石的“三传一反”原理自上世纪50-60年代起就广泛应用于各个学科和工业上。从根本上来说，“三传一反”理论就是宏观层面上的一种说法，其科学意义仅仅知识宏观层面的数学以及屋里归纳。而在实际研究过程中，化学工程在实际研究过程中往往会涉及到较大的时间以及空间维度，时间上涉及分子化学键振动涵盖的纳秒一直到跨度达到数天的整个工业环节，空间上涉及分子纳米单位一直到超过数千米的工厂范围。假设要针对某个具体尺度开展相应控制行为，那么就必须通过另外一个尺度来定位可供操作的办法。时至今日，以往的化学工业体系在过程控制以及放大与优化上还存在诸多问题，而化学工程呈现出往生物学、医药学以及纳米、材料学与环境科学等领域扩充的趋势。所以，从分子层面一直到宏观层面的多尺度联合将会逐步在化学工程领域得到广泛应用。以往的“三传一反”理论就不得不面对更加严峻的考验，这就要求从一个更加新颖的层面来重新认识化学工程领域当中的现象及其规律。与此同时，不均一时空以及多尺度结构仅仅作为化学工程领域内一系列现象当中比较重要的特性之一，而这一点也逐步走进研究者的视野2。李静海等3第一个创建了能量最小多尺度体系，同时在这个基础上成功的预测出快速流化床系统内部某个部分稀密两相的不均匀分布以及密相的聚团大小还有气固流型。这项研究有一个很重要的意义就是表明了用多尺度方法是可行的，也预示对化学工程进行量化也是可行的，并且这种多尺度结构研究将在未来化学工业有一定的发展势头。2.1 多尺度法的基本思路在化学领域中可以将其总结为四中不同阶段：流动、传递、分相和反应；6种尺度：分子、纳微米、单位（颗粒物、液滴、气泡）、聚集体、仪器、生产车间。多尺度法可归纳为：（1）把复杂的总过程分解，从而形成许多尺度各异的子过程；（2）采用不同科学手段对不同尺度子过程进行分析；（3）在研究的过程中，着重于分析不同尺度子过程及其相互关系；（4）针对多尺度形态结构的产生原因利用化学物理方法进行研究；（5）将对不同尺度子过程的分析研究结果进行归纳总结，分析探讨影响总过程的各种因素，从而使总过程中的问题得到解决。可是这些步骤的实施存在很多难题，基本可以分为三类问题：（1）对不同尺度子过程之间的相互联系采用何种研究方法；（2）多尺度多过程进行综合研究时应按照何种原则采用何种方法；（3）如何选择具有代表性的特征量来描述子过程，如尺度空间、尺度等。 通过对不同尺度子过程的研究，我们发现，虽然各种子过程存在差异，但研究过程同样存在共性，如下图所示：图 2 多尺度法研究原理以上多尺度过程和方法可简述如下：(l)分尺度简化。在多尺度结构中会出现一种叫做祸合的现象，主要原因是由于多尺度结构中会出现各种复杂的过程，而且这些过程的发生又在各种不同尺度下。拆分的好处在于子过程中的尺度内部结构相较于总尺度过程在分析研究时要简单的多，所以将一个大系统按照尺度的不同拆分为许多结构简单的组合部分而便于研究。因此在进行总过程的研究时，将总过程进行拆分，从而形成一个个便于分析研究的小过程，也就是说实现了解耦。举个例子来说，聚式流态化中的稀、密两相状态，这种结构是一类较为经典的多尺度有序结构，在这两种结构（稀相和密相）中都有颗粒流体之间的相互作用，但是他们相互作用的情况在这两种相中却是完全不相同的。进一步来讲，就是在密相中，尺度作用是由颗粒进行控制的，而在稀相中这种作用是受流体控制的。两相间还存在聚团和希相两种作用，这两种作用本质上是聚团尺度之间的作用，这种作用是受两相间的相互作用来控制的，两相结构与边界之间发生了设备尺度(宏尺度)作用。利用分解的方法将三尺度进行拆分，使得原本高度不一的结构在分解之后可近似认为是高度相同的稀相、密相以及相互作用，这样通过间接研究这三种相的方法来研究总体就变得相对容易。(2)子过程分析。多种过程之间的祸合是大系统的主要特征之一，在进行分析时，如果首先从总过程着手研究，对认识其内在规律的难度较大，不利于研究的进行，因此只有通过将其拆分为各个子过程，然后将各子过程的研究结果进行综合分析来分析总过程的运行状况。采用这种方法使得研究过程得到了简化，特征子过程发生在各个尺度和不同尺度间的祸合中，采用拆分的方法来处理是基础手段。比如：颗粒尺度作用在流体状态时主要分为悬浮和运送，聚团尺度的发挥造成了运动的无规律性和能量大量消耗。(3)多尺度综合。在以上两个过程都完成以后，为了描述系统的最终模型结构，还需对不同尺度下不同子过程之间的关系进行分析研究，结合已有的条件来进行。这一步骤是三个步骤中最为关键的，因此需要首先确定在进行不同作用间相互耦合时应用何种规则和要求。多值性问题是会常常出现在多尺度系统中的，所以找到确保系统能够稳定运行的前提条件是很重要的。例如：单位质量颗粒耗散能量最大或悬浮输送能耗最小是气固流态化中多尺度作用的祸合原则。一方面微尺度作用中的悬浮和输送过程要求悬浮输送能耗最小，另一方面为维持这一多尺度结构，必然伴有耗散能量最大。2.2 多尺度模拟理论解决较为复杂的工程与材料问题时常常采用的重要办法就是多尺度模拟。这种方法结合了时间与空间上跨层次那个与跨尺度的特点，并且把尺度祸合起来，这样将计算与模拟的效率也提升，从而提炼出有价值的微观数据。多尺度模拟的目标就是抓住不同时空条件下材料或者系统的物理响应特征，并要对材料或系统的性能或者使用寿命进行预测，掌握较小尺度的结构与性能对材料或者系统宏观行为的影响。按照连接尺度的范围大小，多尺度模拟主要包括纳观，微观，细观和宏观等主要尺度的模拟。通常，分子动力学会在微观尺度上进行采用，而在细管尺度上一般采取量子力学理论，当这两个尺度同时采用的时候则采取连续介质力学的防范，不过，与宏观办法相比，细观尺度有一定的不确定性，因而应该与统计学的研究办法结合同时进行计算分析。对于多尺度的建模模拟办法一般有两种方法：第一种是在低尺度层次去进行模拟建模，通过研究，找到内部的法则以及规律，接着，在高尺度的模型中运用找到的规律与法则，该过程是把小尺度上升到大尺度的递增过程。采用这种策略的方法一般称作信息传递的多尺度方法(information -passing multi-scale methods )或递阶的多尺度方法(hierarchical multi-scale methods ) 。还有一种方法是在不同的尺度上同时进行模拟建模，根据不同尺度将区域进行划分而后采用尺度定律进行控制，这些不同的区域能可以灵活选择重叠或者是不重叠，也可以交接连接使用，在这种情况下，处理不同区域间的交接是一个需要重点解决的问题，这种办法通常被称为一致性的多尺度办法。2.3多尺度模拟与计算2.3.1 多尺度问题与方法多尺度问题表现为: 给定某个模型的宏观表述，倘若在局部区域的表述不起作用，则应该相应地调整采用微观非线性这种低尺度的表述来代替。模型的微观特性既受制于宏观上的作用因素，又可能显著影响宏观性能。但微观结构，性能与状态何时、以怎样的途径去影响宏观性能并不清楚，这 些问题对材料与系统等的设计至关重要，同时，也是这类课题研究的热点。假定一个给定系统的微观行为可以使用微观模型变量u 表示，系统的宏观行为用宏观模型变量U 表示，那么宏观模型变量U 与微观模型变量u 可以通过压缩乘子Q 或者重构算子R 联系起来: ( 4 ) ( 5 )如果在整个系统中使用微观模型，建模太复杂，计算量太大。与此同时，相较于系统的宏观性特性，微观性则显得重要性不那么高。例如，对于输送管道系统而言，重要的是管道每天的输送能力如何，至于管道中的流体分子在管道中如何运动并不重要。因而，需要将微观模型作为宏观模型的补充,或者划分出局部区域建立微观模型，使之为获取系统宏观行为提供必要的信息，从而实现系统宏观行为的有效地多尺度描述，显然要比在整体上使用微观模型更加高效，也更容易实现。2.3.2 多尺度模拟与计算方法采取的计算办法与多尺度模型应依具体问题而定。就多尺度建模和计算而言，可以将多尺度问题按其自身特性进行分类，根据每一类情况的特性，又发展出不同的建模与计算方法。本文就一般常用的两种多尺度计算与模拟问题进行介绍。一类是含有分散的孤立缺陷问题，另一类是基于微观模型推测宏观行为问题。2.3.2.1 含有分散的孤立缺陷问题对于含有分散孤立缺陷这一类问题(固体中的微裂纹，位错，振荡，流体中的接触线(contact line)以及生物多聚体中的酶(促) 反应等)，其特点是: 建立模型时，需要在出现缺陷或者奇异性的局部区域建立微观模型，而远离这些缺陷或者奇异性的其它区域仍然可以使用宏观模型。对于常见的一些缺陷，如单晶体中位错与微裂纹等，可以采用非局部准连续办法进行模拟。在此情况下，微观模型采用分子动力学模型，宏观模型采用非线弹性模型，从而实现纳观一连续介质力学和微观一连续介质力学间的连接。自适应的网格细化方法使得准连续体方法可以识别出缺陷的局部区域，并能够局部细化到原子尺度，以期求解包含缺陷的整个局部区域的详细信息。但忽略了原子的振动，因而在低温时静态模拟精度较好。系统总体的能量计算公式如下： (6)上式中为代表性原子的能量，为代表性原子的总数目，为由代表性原子代替的原子数目，其中E为系统总能量。 由Abraham等人提出的MAAD (MacroAiomistic ab-initio Dynamics) 方法，对随后的许多其它多尺度模拟方法的发展有很大的促进作用。该办法首先用在硅的裂纹扩展仿真研究当中，采用最基础的TIght-Binding模型对裂纹的断裂键能进行模拟。在围绕裂纹尖端区域使用分子动力学模拟位错现象，远离裂纹尖端区域使用有限元模型，给动力学施加边界约束，从而获取宏观上的形变。这种方法实现了远离裂纹区域有限元模拟，裂纹尖端区域的分子动力学模拟和裂纹尖端处的Tight-Binding 模型之间的无缝祸合。这样的仿真办法结合了连续介质力学、分子动力学以及量子力学三种方式进行整合仿真。这样得到的系统的总体的能量哈密顿计算公式如下： (7)上式中，C是刚度，为应变量， ，是材料密度，则是第i个原子的速度，V为MD区的总体势能，为原子间相互排斥项和相互吸引项的共同作用。过渡区的能量(， ) 将过渡区两边的能量函数进行线性结合来选择，系数的选择依据每个过渡区的函数选择的贡献能量。各区域相互传递能量的多少可以忽略，可以采用同时间步的方式计算各自能量函数。把单元间原子级的范围降低能够增加计算的准确性，这样计算成本与计算量也会增加。在裂纹扩展的扩展仿真模拟上这种办法也运用得比较广泛。CGMD (Coarse Grained Molecular Dynamics)方法由Rud和Broughton等人提出。其基本思想是采用MD方法计算需要详细描述的关键区域，而引入有限元方法粗化计算远离关键部分的区域，即挑选一定密度的原子作为节点，将其他原子只作简化计算，以此降低计算消耗。远离关键区域的粗化能量简化表达式为: (8)其中表示从系统中被粗化掉的原子的能量，当节点数目接近原子数，即接近关键区域时，粗化的原子能量可以被忽略，对MD计算起主导作用。为动能表达式，为势能表达式。质量矩阵和刚度矩阵可以借助有限元方法的形状函数求得。本质上CGMD方法与QC方法的思想类似，但是采用了更加有效的计算公式所以使计算量被减小了，而且解决了温度问题的也是其优点之一。然而粗化能量公式只能计算一些能量公式可以简化的特殊情形，原形实际上并不易计算，这一点限制了它的使用范围。CGMC (Coarse-grained Monte Carlo ) 法是一种半解析半数值的多尺度方法，使用长程相互作用系统的粗粒化蒙特卡罗模型。通过适度的平均化微观模型中的参数，由解析方法得到粗粒化模型的相互作用参数。在粗粒化水平强制细致平衡，微观模型和粗粒化模型处在同一个细观尺度保证了系统的自洽性。CGMC法在复合材料和流体等方面有广泛的应用。DSMC (Direct Simulation Monte Carlo ) 法由Bell，Berger，Colella 等人近年来发展起来的。作为自适应网格细化法的推广，在振荡局部区域采用较低尺度的动力学模型，其它区域采用连续介质力学模型，并利用离散仿真蒙特卡罗法求解动力学模型。通过元胞边界变换通量将连续介质力学模型和离散仿真的蒙特卡罗法祸合起来。用于流体的连续体一分子动力学模型在流体，乃至微射流的模拟中，研究者们对分子动力学与连续介质力学祸合模型进行了大量的尝试。其基本思想是异构的域分解: 将求解域划分为重叠的区域，在这样的区域上进行连续体和分子动力学模型的求解。要求连续体和分子动力学模型的结果在重叠区域相匹配，以保证整个计算的连续性。用于固体的连续体一分子动力学模型将连续体分子动力学模型祸合起来以模拟固体的动态响应特性。存在的主要问题是连续体模型一分子动力学模型结合处的匹配条件。Cai 等人建议在处理线性问题时使用真实的边界条件，但是，这样做的计算成本很高。E .和Huang 等人将传统上用于波的连续体计算的吸收边界条件推广到了分子动力学的离散模型中。这可以算是真实边界条件和吸收边界条件在连续体模型中的直接简单推广的很好的折中。这种方法在处理低温，系统无热能输入的情况下相当有效，在模拟裂纹扩展中也相当奏效。2.3.2.2 基于微观模型的本构模拟问题这类问题的特点是: 无法使用宏观模型对系统进行有效地模拟，应该间接地通过较低层次的微观模型推断其宏观行为。在分子动力学的模拟中，原子势能不明确，或者对固体和流体进行连续介质模拟时，没有明确的本构关系，这种情况就会出现。通过直接利用原子势能，而非经验的贮能函数，使准连续体法同样能够用来模拟材料的宏观变形。由于一小簇原子的能量可以用来表示整个宏观元胞的平均能量，计算就不需要深入到大的空洞里的原子，比计算系统全部的原子更加高效。系统总能量计算公式为: (9)其中E 为系统总能量，为示单元数，为单元体积，即任意单元中的变形梯度不变。人工压缩法用于不可压缩流体的模拟计算。在这种情况下，利用可压缩流体的方程模拟计算可压缩流体，却不使用不可压缩流体的方程。尽管不可压缩流体的方程Navier-Stokes方程是可知的，但其数值求解不容易实现。通过合适地选择马赫数(Machnumber)，在计算精度和计算成本之间做出一个可以接受的折中。气体动力学方法。尽管发展这一方法的初衷并不是为了多尺度建模，但却为多尺度模拟提供了最好的实例。这一方法的目的是直接使用动力学方程而不是气体动力学的欧拉方程来计算气体的宏观行为。采用重构一演化一平均化的方法，这一思想成为异构多尺度方法(HMM)的起点。HMM法由WeinanE 和Engquist B等人提出。采用自上而下的建模策略: 先建立一个非完整的宏观模型，再将微观模型作为宏观模型的补充部分，实现不同尺度间的祸合。HMM法有两个主要组成部分: 宏观求解器和从微观模型估计数值信息丢失的过程。其中宏观求解器是基于已有的知识,即它的设计或选择应该尽可能多的考虑人们对宏观模型的己有认识。而数值评估这一步的关键是能够找到一个合适的有约束的微观求解器。这种约束的必要性表现在它能够确保微观模型的结果与系统宏观模型的结果相符合。而这是HMM法最具有挑战性的一步，并与所要解决的问题紧密相关宏观模型与微观模型整个祸合过程通过如下方式来完成: 系统的宏观模型为微观求解器提供约束条件，微观模型向宏观模型求解提供必要数据(通量,受力情况和刚度矩阵等)，这种基于数据的祸合策略给多尺度建模提供了非常大的灵活性，在多个领域中的都有着广泛的应用。需要指出的是，这种方法与前面提及的其它多尺度建模与计算不同，它不是针对某个具体的多尺度建模与计算问题提出的，而是试图建立一个设计多尺度建模与计算的通用的指导框架，最大限度地利用人们对系统宏观行为己有的认识，实现高效地多尺度建模与计算。因此，可以认为HMM旨在指导设计高效的多尺度建模与计算方法。 3 多尺度模拟与计算机关系通过计算机仿真模拟建立一个简单的模型，可以找到符合预期性能的最佳配比以及最合理的工艺流程，这样可以节省大量的人力、物力和财力。凭借巨大的存储空间、高速运算能力和逻辑判断能力，计算机在多尺度模拟中所产生的作用是至关重要的，提供巨大的帮助，因此计算机模拟在研究中扮演着构筑理论方法和实验方法桥梁的角色。计算机会通过大量运算对设计方案进行优化和综合分析，最终确定最合理的设计方案。通过这条途径会使设计质量显著提高，设计周期得以大大缩短，从而产生的效益就事开发成本的显著降低。计算机会激发出人类创造性的构思方案，通过对存储资料进行筛选以及方案权衡比较，找到最佳方案。此外，一些现有物理和化学实验无法实现的极端条件，比如高温等离子体、核反应堆、高温高压、高辐射条件，可以利用计算机实现对等进行模拟研究。由于高分子组成多样，又具有长链结构，其物理和化学性质有特殊性，在时间和空间表现为多尺度性。例如高分子嵌段共聚物经过微相分离可形成尺寸为几百个埃的特征相区，又比如高分子主链单键键长虽然为埃量级，但它的单链回转半径可达到单键键长的1O0倍，所以显然不可能单单的从原子级别的模拟来研究高分子的性质。而动力学在时间尺度上跨度更大，比如虽然键振动一般为1Os1 5S量级，而高分子结晶和增长过程所需反应时间则常常超过了1S，高分子共混物的相分离甚至需要更长的时间完成。高分子体系的模拟研究显然无法通过单一的计算机模拟方法来完成。传统的计算机模拟方法主要包括蒙特卡罗方法，分子动力学方法和布朗动力学方法等等。其中第一种通过对构型空间进行抽样从而来计算系综平均的基于格子模型的蒙特卡罗方法是最常用的一种；第二种分子动力学则是根据经典力学原理，使用力场的研究方法，对粒子运动的动力学方程求解，模拟系统随着时间的推移演化。最后一种布朗动力学中通过Langein方程描述粒子的运动方式，布朗动力学的主要特点在于考虑了体系中粒子的热运动带来的随机力。分子动力学是一种非常重要的计算机模拟方法，它主要研究对象是原子和分子，是对多体问题进行求解的方法与技术，被广泛的应用于多粒子体系的研究中。蒙特卡罗方法是一种统计学上的概率性方法，最常用的是重要性抽样来研究相空间中马科夫链的随机行为。如果我们对于某个体系(如聚合物溶液)中的溶剂分子的具体运动不感兴趣，耗散力和随机力就可以取替溶剂对聚合物的作用，这样分子动力学中的牛顿运动方程就会转变为郎之万方程，形成了布朗动力学的模拟方法。以上这三种方法的模拟尺度在时间上为纳秒范围，空间上不足一百个埃，因此可被归为微观尺度的计算机模拟方法。另一方面介观尺度的格子玻尔兹曼，动态密度泛函，耗散粒子动力学以及场论聚合物模拟等方法和理论得以大发展，因为本质上，尺度与高分子材料的物理特性直接相关而且往往体现在微秒甚至毫秒这种级别。随着人们对材料宏观性能需求的增加，以有限元为代表的宏观模拟方法应运而生。 至今，人们已经实现了从宏观细粒化到介 i，进而到微观以及经粗粒化实现反方向这样一个贯穿的模拟过程。介观的耗散粒子动力学模拟方法忽略分子或者链上得某些可以忽略的信息，经粗粒化微观尺度上得到的几个分子或者高分子链就构成了耗散粒子动力学中的粒子。耗散粒子动力学在研究复杂流体方面具有很突出的优势，它作为介观尺度的计算机模拟技术已经应用到了研究方面，比如化工、生物和医疗等诸多领域。从耗散粒子动力学的本质上来说，它属于离散粒子动力学方法，这就导致耗散粒子动力学处理的体系不能太大，从而影响到它展现真实流体的物理特性。离散颗粒模型（DPM）在颗粒流体系统的模拟中已经得到了广泛的应用，但离散颗粒模型对大规模的系统进行模拟时计算量巨大，这种巨大的计算成本严重制约了它的应用和推广。然而技术上不断发展突破的的图形处理器（GPU）以其强大的计算能力和良好的并行性为离散颗粒模型的发展应用提供了另外一种可行性。本文利用CPU 计算颗粒相和流体相的运动方程建立的CPU-GPU 耦合的计算模式，基于这种模式模拟了颗粒流体系统随着设备结构的尺寸的变化而变化，从节涌、湍动流态化向快速流态化过渡。在此基础上，本文进一步实现了CPU-GPU 耦合的多尺度并行计算模式，模拟研究了包含有500000 个颗粒的气固流态化系统的非均匀结构的动态演化过程。模拟结果复现了上稀下浓、中间稀、边壁浓的稀密相共存的环核结构的形成过程，展示了CPU-GPU 耦合的多尺度并行计算模式在颗粒流体系统模拟中的应用前景。由于颗粒轨道模型能够在微观尺度上使我们获得单个颗粒的具体信息，这个功能对于研究颗粒流体系统中复杂的流动行为提供有效而深远的帮助，同时在实际工程操作与应用中表现出巨大的潜力。但仍然存在的一个现实问题是颗粒轨道模型的应用依然需要大量的计算时间，这是因为实际工业装置的规模都比较大而且具有复杂的边界，同时实际的工业过程通常伴有复杂的流动和传热以及其他化学物理反应。近年来，计算机图形处理的速度和图形处理的质量通过图形处理器(GPU)的迅速发展被极大地提高了。同时，图形处理器（GPU）在绘制流水线是呈现了的很高的速度和良好稳定的并行性，甚至近几年来逐渐发展起来了GPU的可编程功能，不仅为图形处理，甚至包括数字图像处理和通用并行计算都因此得到了良好的运行平台，使得人们越来越关注如何使用GPU来并行加速图像处理算法和传统的通用计算。至于时间成本最大的颗粒轨道模型中的颗粒相间作用，也利用 GPU 进程间快速便捷的通信，减少传统并行机进程间通信和管理的损耗；同时因为粒子规模与GPU 的执行时间是次线性的关系，利用GPU 的并行特性和高速浮点数计算不仅提高了算法的运算速度，还可以较大幅度地增加计算的规模。针对传统并行颗粒轨道模型在实际应用中的不足，本文结合图形处理器(GPU)的高速并行性，提出了一种基于GPU 加速的多尺度并行颗粒轨道模型，在解决了计算的粒子规模问题的同时尽量减少了需要耗费的计算时长，从而使得颗粒轨道模型的运算速度得到提高。目前，解决这一问题有两个办法，其中之一是可以通过提高CPU的运算能力，以及借助GPU运算来实现计算机性能的提高。其二就是采用高性能并行计算技术。通过采用并行计算技术可以将复杂系统的庞大计算任务分给多个CPU共同完成，扩大了模拟体系使得模拟结果结果也会更加准确，从而展现真实体系的物理特性。这种办法会显著提高计算效率因为参与通信的粒子数量是非常少的，如果与被分配到各个微处理器上的粒子数量相比的话。开展高性能并行计算程序、提高计算机水平、借助GPU强大的处理能力和良好稳定的并行特性、发展GPU通用运算研究，势必会对在较大的时空尺度上的复杂流体的研究产生重要的作用。随着模拟方法的不断完善、并行技术的不断发展和性能的不断提高，计算机模拟在人们未来研究过程中将会起着重要的纽带的作用。3.1 CPU简介所谓中央处理器（CPU），其主要是指超大规模的集成电路，通常情况下，其是计算机的运算核心与控制核心。其主要功能有两个方面，一是对计算机指令予以解释，二是对计算机数量进行处理。CPU主要由以下几个部分组成，一是高速缓冲存储器，二是控制及状态的总线，三是运算器，四是实现他们之间联系的数据。CPU与输入/输出及内部存储器合成为电子计算机三大核心部件。对于控制部件而言，主要完成译码指令的工作，译码之后还需发送相应的执行命令。其结构有以下两种：1、以微存储为核心的微程序控制方式；2、是以逻辑硬布线结构为主的控制方式。一条微码则对应一条微指令操作，该微码存储在微存储当中；一条指令是由许多个序列不同的微码构成的，这样的微码序列由形成一个完整的微程序。这些指令在CPU当中译码成功后，被转换成一定时序的控制指令，根据微码的序列，依次发送，其执行的节拍则由微周期决定。也可称为指令的执行。CPU在高速缓冲储存器或者一般储存器中读取操作指令，而后将其传送至指令寄存器，在通过译码后，这些操作码会进一步分成有顺序的微操作，最后成为可执行的命令，这样便是执行指令的完整流程。计算机执行或者操作的基本命令就称为指令，它由若干的自介构成，其中可