应用物理学毕业论文.doc

西华大学毕业论文摘 要格子Boltzmann方法从流体的分子运动论出发，采用统计平均的方法来研究流体的宏观运动和性质，同时它还是能与高性能计算机相匹配的新型流体力学计算方法。本文用Fortran语言，MPI并行算法编程模拟牛顿血液流的流场分布，使用的计算方法为三维格子Boltzmann方法，模拟人体血液系统中心脏动脉弓中血液的流场分布，将计算结果用图形的方式绘出，分析动脉弓横向截面中血液压强的变化规律。以此分析心血管疾病研究动脉粥样硬化与血液动力学因子的关系，为临床医学提供参考。关键词：格子Boltzmann方法 横截面 动脉弓 压强AbstractThe lattice Boltzmann method base on molecular motion of the fluid, use the statistical method to study the nature of fluid .At the sometime it is a calculation method which able to model the fluid mechanics and high performance computer matching. This article using Fortran language, MPI parallel programming algorithm to simulate the distribution of the Newton flow fluid by the 3D lattice Boltzmann method. We simulate the flow of the human blood in an aotic arch, the calculation results are drawn by the tecplot. Through those picture, we analyis the change of blood pressure in arterial arch tranverse section under the Reynolds. Using the lattice Boltzmann method to analyze the relation between atherosclerosis and blood hemodynamic factor can provide the reference for the clinical medicine.Keywords: the Lattice Boltzmann method; transversal surface; aotic arch; pressure. 引 言随着计算机技术的高速发展和信息社会的到来，对各种复杂现象的仿真模拟渐渐成为了继实验和理论计算后又一重要的科学研究工具，越来越受到人们的重视。而以解决大规模科学与工程问题为背景的高性能计算更是其中的佼佼者。近年来，人们逐渐认识到，高性能计算关系到国家在现代科学和高科技研究与开发中竞争能力和领先地位的问题，因而受到当今世界各国的普遍重视和关注。1999年2月，美国公布了TT2计划(21世纪美国信息技术)，将高性能计算和先进的科学工程计算列为重点支持方向，并指出，人类已经进入科学计算作为科学和工程研究的手段，与实验和理论同等重要的时代，21世纪的复杂问题将通过科学模拟来解决。而其中流体力学的仿真模拟是大规模科学和科学计算的一个重要领域，尤其是随着流体现象复杂性的增加，如何使用高性能计算解决流体问题，成为大家亟需解决的问题。流体是由大量分子构成，每个分子都在不停息地做无规则的热运动，并且通过频繁的碰撞交换动量和能量，而流体的宏观运动和其它性质是流体分子微观运动的平均结果，与此相对应，研究流体运动有宏观和微观层次上的两种途径：宏观上流体的运动以连续介质为出发点，可以用一系列微分方程描述（如Euler方程、Navier-Stockes方程等，简称N-S方程），这些方程大多是非线性的，要想得到精确的解析解是不容易的。故人们常用的方法是在时空域中将N-S方程离散化，根据对方程中流变量所采用的离散化方法不一样，可将这些不同方法分别称为有限差分法、有限元法、有限体积法、谱分析法等，然后对离散的方程求数值近似解。目前人们开发出来的可以做以上工作的专业流体力学模拟软件有FLEUNT、FIDAP、CFX、PHOENIX和STAR-CD等， 其中计算流体力学（CFD）应用最为广泛，已经应用于航天航空、大型能源装置（核电站）、新型交通工具、环境保护等众多技术部门和领域，成为现代计算科学不可或缺的工具。但是随着运用的进一步推进，它的局限性也逐渐暴露出来，例如，由于计算条件的限制和缺乏有效的计算模型，CFD对航天器热问题的计算误差可以高达100%。所以，科技的发展对计算流体力学提出了更加严峻的要求。在这种历史环境下，从微观角度来研究流体力学的新方法进入了人们的视线，这种方法就是格子Boltzmann方法。 格子Boltzmann方法从流体的分子运动论出发，采用统计平均的方法来研究流体的宏观运动和性质，同时它还是能与高性能计算机相匹配的新型流体力学计算方法。据悉，我国已经研制出曙光2000和神州3000这样的大型机，并且每秒千万亿次的并行机已经面试，并很快推广起来，但是无论国内还是国际上，由于模拟计算的算法的限制，并行机的运用率并不高，新的并行机需要新的科学计算方法。为了有效地利用并行机，我们需要的算法是要把问题分解成在不同计算机上进行处理的多个部分，然后重新组装在一起，这些正是易并行的格子Boltzmann方法1所能完成的。故自问世以来，格子Boltzmann方法引起了物理学家、数学家以及计算机科学家们广泛的兴趣。格子Boltzmann方法不但是晶格气体自动化方法的发展，也是连续Boltzmann方程的一种特殊有限差分格式。格子Boltzmann方法是以有限速度集上的粒子分布函数为基础，实质就是将空间和时间均离散化，空间分为均匀的网格，每个网格单位为x，时间分为均等步长t。粒子只能位于节点上，在节点上与从相邻节点运动来的粒子相遇，并发生碰撞。碰撞后，粒子以碰撞后的速度运动到其它节点，且通过碰撞算子对粒子分布函数进行演化。在演化中根据节点上所必须满足的质量、动量和能量守恒，宏观量如质量，压强等均可通过粒子分布函数求得碰撞算子可以选择多种形式，但是单弛豫时间近似的碰撞因子最为简单，故人们趋向于使用这种形式的碰撞因子。单弛豫时间近似的碰撞因子虽然是非线性的，但是通过Chapmann-Enskog方法对格子Boltzmann方程中的时间导数、空间导数和分布函数进行多尺度近似展开，可以在时间步长的二级精度下恢复宏观的Navier-Stokes方程，实现了从微观Boltzmann方程到宏观Navier-Stokes方程的转化，从而证明了可以使用格子Boltzmann方程来模拟流体的运动，且流体宏观的压强和切应力均可用粒子分布函数直接求出，避免了宏观传统计算方法中求解压力需要求解Poisson方程和切应力需要求速度的偏导数的麻烦，这就是为什么说格子Boltzmann算法简单的缘故。1研究现状与传统的CFD相比，格子Boltzmann方法具有算法简单、计算效率高，边界条件易处理，容易实现并行等优势，吸引了众多学者对其进行研究。在二十年的发展过程中，格子Boltzmann方法被越来越广泛地运用于各种复杂流体的计算1，诸如粘弹性介质2，粒子沉降及悬浮粒子在流体中的运动3，多相流4，化学反应扩散5、非牛顿流体6等。热格子Boltzmann方法还能模拟 Benar 花纹7,8，格子Boltzmann方法在研究湍流方面9也取得了成功。目前，格子Boltzmann方法的研究仍在不断的发展中，每年都有新的模型产生，其应用领域也在不断扩大。与此同时，格子Boltzmann方法在血液动力学中的应用也取得了长足的发展。血流动力学以血液与血管的流动和变形为研究对象, 将力学的理论和方法与生物学、医学的原理和方法有机地结合起来，研究血液和血管的力学特性，分析血液在心血管系统中的流动规律，探讨某些血液循环系统疾病对血液流动特征的可能影响。目前，国内外很多单位都针对将格子Boltzmann方法应用于血液动力学做出了很多尝试，例如 Mittal 等人就用 large-eddy 模拟方法研究脉动流在有狭窄的动脉血管中的流动10， Cristiano 等人和 Sun 等人用 LBM 研究了毛细血管中白细胞对红细胞的运动的影响11,12， Pozrikids 用 Green 函数法研究了红细胞在剪切流中的运动13,14， Secomb 等人用润滑理论研究了红细胞在毛细血管中的运动15,16。国内，中国科学院上海应用物理研究所的方海平教授等人为了研究血液流，先后提出了基于晶格Boltzmann方法的数值模拟可扩张的管道内的静态和脉动流、单个小颗粒在流体中的运动等的计算模型，以及在红细胞变形等情况上取得了很大的成功。应用三维格子Boltzmann方法对动脉弓模型的Casson血液流场进行模拟。管道弯曲是动脉粥样硬化局部性的另一种代表性模型，。在不同雷诺数情况下进行了模拟，给出了弯道中血液流流场的分布，二次流的发展情况，确定了在动脉弓的弯曲管道中可能发生动脉粥样硬化的部位。我们课题组主要针对运用格子Boltzmann方法研究动脉粥样硬化与血液动力学因子的关系，为临床医学提供参考。动脉粥样硬化（atherosclerosis）是一组动脉硬化的血管病中常见的最重要的一种，其特点是受累动脉病变从内膜开始。一般先有脂质和复合糖类积聚、出血及血栓形成，纤维组织增生及钙质沉着，并有动脉中层的逐渐蜕变和钙化，病变常累及弹性及大中等肌性动脉，一旦发展到足以阻塞动脉腔，则该动脉所供应的组织或器官将缺血或坏死。由于在动脉内膜积聚的脂质外观呈黄色粥样，因此称为动脉粥样硬化。20世纪70年代以来,随着新技术的应用,促进了从血液动力学角度来探讨动脉粥样硬化的发生发展和定位的机理，并已经知道，动脉粥样硬化和其它动脉损伤一样，并不是随机地出现于系统脉管中，某些动脉几何形状对促成这种疾病有着相当高的几率，动脉粥样硬化的流体力学理论就是在这种基础上建立起来的。格子Boltzmann方法作为模拟流体运动的一种数值计算方法，它不但是格子气自动机5（Lattice Gas Automata，简称LGA）模型的发展，也是连续Boltzmann方程的一种特殊的有限差分格式，实现了从微观Boltzmann方程到宏观Navier-Stokes方程的过渡。本章从格子气自动机和连续Boltzmann方程两个方面介绍了格子Boltzmann方法的基本原理和基本思想，并利用多尺度展开将格子Boltzmann方程恢复成宏观流体力学的Navior- Stokes方程，证实了格子Boltzmann方法模拟流体的可行性；最后给出了格子Boltzmann方法常用的速度离散模型，并给出了该方法的计算框架和算法。2.压强边界和曲线边界条件除了固边界，我们还要处理入口 和出口的边界情况。在本文中，我们对于入口和出口主要处理了流体在有压力差的驱动下的流动（压强边界条件）和边界有滑移时的流动（速度边界条件）情况。Zou 和 He 提出，在反向弹回边界条件中，在满足质量守恒的前提下，修改各个边界格点上的分布函数，使之满足一定的流体压力或速度。以D2Q9模型为例来说明，如图2-1对于出口和入口的任意一点，跃迁后根据密度和速度的定义式可得： (2.1)下面的压强和速度边界条件都是基于方程组(2.1)来介绍的。2.1压强边界条件压强边界条件就是给定入口和出口的压强下求未知的分布函数，我们知道给定压强也就是给定密度。对于入口处的节点，如图2-1在跃迁后，已知，未知。令，而未知。由方程组(2.1)可以把解出来： (2.2)图2-1 压强边界条件2.2曲线边界条件如图2-2所示，一条曲线边界存在于格子单位长度为的格点之间，如果一对格子沿着跃迁方向有一个以上的连接与曲线相交，那么流体内部的格点用来表示，边界外部的用来表示，与曲线边界相交的实心点用来表示，令其速度为（对于非滑移，）。所谓边界条件，就是要把从边界格点流到流点的分布函数确定出来。由图中我们可以定义 图2-2 曲线边界条件示意图3.物理模型在血液的格子Boltzmann模型中，我们使用Fortran语言编程模拟。在研究早期，为了计算和建模的方便，我们认为在动脉血管中，血液的有形元素相对于动脉管径的大小可以忽略不计，故将血液近似为牛顿流体。在此近似下，我们课题组前期进行了血液流在血管中流动的模拟。我们先来看看所采用的格子Boltzmann方法的算法。格子Boltzmann方法中粒子在每个节点上的演化分为以下两个步骤：碰撞过程：在节点上的粒子与从相邻节点运动来的粒子发生碰撞，根据质量、动量和能量守恒规则改变粒子的速度: (3.1)迁移过程：碰撞后的粒子从节点在一个时间步长内以恒定的速度 沿着方向运动到相邻的节点，再与其上的粒子发生碰撞： (3.2)算法3-1根据血管形状，写边界函数bound初始化分布函数：while(满足终止条件-流场收敛)do （使用迭代方法进行循环）for each node , do碰撞: endfor each node , do迁移： endfor each node , do计算宏观物理量： endend输出计算结果 以上两步交替循环直到流场达到收敛为止。编程中流场收敛的判据我们取为： (3.3)表示第n次循环后的粒子分布函数，上式求和在整个流场中进行。格子Boltzmann算法。如算法3-1所示，一般地，我们认为在初始时刻流场是静止的，故初始的粒子分布函数取速度为0的平衡态分布函数。程序的流程可以用图3-1来表示。辅助步骤物理建模，确定计算区域，初始条件、边界条件等计算区域网格划分、确定节点、边界点全场初始化，确定各节点上的宏观参量即粒子分布函数模拟步骤跃迁边界处理碰撞计算各节点的宏观物理量，输出计算结果判断是否收敛是否图3-1 格子Boltzmann方法的编程流程4.三维格子Boltzmann模型的MPI并行编程设计 流体力学的计算与仿真问题一直是高性能并行计算研究与应用的最重要领域之一，但是由于传统的流体力学计算基于一组耦合在一起的非线性宏观偏微分方程，因此在运用于并行时遇到了很大的困难，并行算法在传统流体力学计算方法中基本上是“英雄无用武之地”。直到元胞自动机（LGA）理论的出现，诞生了流体仿真模拟的完全并行的格子气自动机，为流体力学的并行计算模拟提供了实现途径。而格子Boltzmann方法继承了LGA方法优点，同时比起LGA来说，碰撞因子更加简单，也不需要进行统计平均，因此具有更好的并行性，更适合进行并行计算。图4-1 MPI并行程序数据处理过程并行程序由一组独立运行的进程构成. 进程间通过相互发送消息来实现数据交换。在使用MPI并行程序时，首先将编程区域分解成多个子区域，分配给不同进程，保证进程间的负载平衡和最小的消息传递通信开销；然后选择合适的通信数据结构，这里我们使用的是Jacobi迭代，其局部性好，可以取得很高的并行性；在MPI程序中，为了并行求解，将参加迭代的数据按其结构进行分割，各块中相邻的元素进行通信，内部则完全独立并行计算，当计算规模增大后，通信的开销相对于计算的比例会降低，这样使用并行计算的优越性就显现出来了。为了简单起见，我们以二维问题来说明一下Jacobi迭代的数据划分和通信形式。假设需要迭代的数据是的二维数组A（M，M），将数组A进行分割，平分成四份，即有四个进程同时进行并行计算，每个进程使用的数据为二维数组B（M，N），划分结果如图4-1所示。每个进程上加载的数据分别是：进程0：A（M，1：N）；进程1：A（M，N+1：2N）；进程2：A（M，2N+1：3N）；进程3：A (M, 3N+1: M)。在迭代中每个进程边界点的计算需要相邻进程的数据支持，因此我们在每个进程数据块的两侧各增加一列数据的空间，用来存放相邻进程传递过来的数据，这样每个进程上加载的数据块大小变成了。数据块划分好后，我们来看一下数据的通信和计算流程：在迭代之前，每个进程都需要向相邻的进程提供数据块，同时接收相邻进程传送过来的数据块，具体传输过程如图4-1所示。我们使用二维格子问题来阐述如何设计基于区域分裂的并行格子Boltzmann算法，假设格子Boltzmann方法使用的网格整体规模为，并假设将流场网格分裂成 H个大小相等的子网格，则，使每个进程（处理器）处理 个网格节点，网格分裂方式如图4-2。在格子Boltzmann方法中，粒子分布函数仅仅在迁移步骤与它的相邻节点发生联系，即 (4.1)图4-2 网格分裂示意图 所以在每个进程中对所有的内节点进行迭代计算，而边界节点上的值，我们需要知道与之相邻的子网格的边界节点的值，为了方便存储相邻子网格传送过来的数据，我们在每个边界层外面再加两层，使每个进程上的网格数变为，这两层网格上的相关数据由相邻的子网格获得，同时该子网格中的边界节点上的数据也要发送给相邻的子网格。因此在每次迭代之前，相邻的进程间需要先互相交换信息，信息交换模式如图4-1所示，然后在每个进程中所有节点进行跃迁。跃迁之后，对边界点执行相应的边界条件，对所有节点执行碰撞步骤。如果不满足结束条件，接着在各进程间再进行数据交换，然后进入下一次迭代，直到流场收敛为止。5.动脉弓的血液流场模拟众所周知, 在弯曲动脉中容易发生动脉粥样硬化等心血管疾病。动脉粥样硬化的这种病灶性是与弯曲动脉中的血液流动特性密切相关的, 即与血液流动的力学因素和动脉血管几何(形状) 因素的相互作用机制有着极其密切的关系。血流速度和压力的改变、流动的分离与再附、二次流的产生、血管壁面切应力的变化，造成了血管的局部生理变化,，导致了动脉病变。弯管动脉内的血液流动问题是生物流体力学的一个重要研究课题。5.1动脉弓模型作为数值模拟的几何图形，我们就将动脉弓简化为曲率半径Rc = 6 mm、环向弯曲为180的弯管，且血管各截面内半径R = 2 mm ，。图5-1显示被离散后的三维动脉弓的对称二维横截面内的网格划分及边界点分布，为动脉弓内某一垂直于轴的横截圆面与水平面的夹角，血液由弯管的右边管道流入，对于胸主动脉弓来说，右边的直管部分，称为升主动脉，其与心脏相连接， 和动脉弓一起成为动脉系统的开端，而相应出口的直管则称为降主动脉。我们课题组之前使用三维圆直管进行了三维格子Boltzmann模型优化的研究(过程类似4.1节)，最后确定了在模拟中选用的格子Boltzmann模型，即D3Q19模型。边界条件方面，管的入口和出口的速度均被认为全面发展，并且让入口和 图5-1三维动脉弓在yz面上的截面（黑点为边界点）出口分别处于对应节点的连线中间，血管壁无滑移；入口、出口以及管壁都采用Mei65 等修改的FH曲线边界条件。模拟中，出口和入口的密度为：，空间步长取为1mm，格子数明显超过单个计算机的计算负载量，故我们将网格沿着z轴平分为8份，进行包含8个进程的并行计算。5.2模拟仿真三维动脉弓，建立模型并分析模型1.运用fortran语言，模拟仿真三维动脉弓，建立模型。并通过MPI并行运算，得出在三维心动脉弓中血液流场的压强分布的数据库。2.运用fortran语言，编译程序，实现对数据库的处理，得到在三维动脉弓从入口出发到出口过程中0度、30度、60度、90度、120度、150度、180度不同圆截面上的压强分布数据。3.运用tecplot软件，对数据进行处理，绘制出在雷诺数R=760下不同圆截面中压强的等高图。如图5-25-5所示，图中上面为外侧壁，下面为内侧壁。图5-65-8为上面为内侧壁，下面为外侧壁。 图5-2 U型血管0度时 图5-3 U型血管30度时 压强的横截面等高图 压强的横截面等高图 图5-4 U型血管60度时 图5-5 U型血管90度时压强的横截面等高图 压强的横截面等高图 图5-6 U型血管120度时 图5-7 U型血管150度时压强的横截面等高图 压强的横截面等高图图5-8 U型血管180度时压强的横截面等高图图5-2为在雷诺数为760时U型血管入口处直管0度时的压强横截面等高图，通过观察可知开始血液进入血管时外侧压强最大，高压部分比较多。而内侧的低压部分比较少。而且从外侧到内侧压强均匀地降低，压强梯度基本恒定，又遇外侧压强大，内侧压强小，圆截面处血液由外侧流向内侧并且沿直线流动。就说明在血液刚刚进入血管时外侧血管壁压强较大，内侧较小。图5-3为在相同雷诺数下U型血管30度时的横截面压强等高图，通过和图5-2比较不难发现高压部分明显减少低压部分明显增多。而且压强从外到内基本均匀降低但是入口处要乱一些。这就说明了在30度时外侧血管壁的压强比刚刚进入血管时有所减小，高压区域减小。图5-4则是60度处进入弯管后的横截面压强等高图，经过比较我们可以很容易的看出高压部分更加的少了，低压部分已经占了很大的比例，这就说明了此时的外侧血管壁压强比30的时候又有所减小。压强梯度变化尤为剧烈不在均匀，此时还是外侧压强大于内侧压强，血液由弯管外侧流向内侧，但是流动轨迹不再是直线，而是曲线，支持血液偏离圆面中心，向圆周处流去，形成二次流。图5-5为90度时横截面的压强等高图，通过观察我们会发现高压比前一幅图有所增加，可是仔细看一下等高图表的最大值，我们会发现不管是外侧血管壁还是内侧血管壁的压强值都有所减少，所以压强也在减少。血液压强梯度进一步变化剧烈，血液流场二次流明显。图5-6是我们取横截面为120度处流出弯管后U型血管中压强的等高图，压强梯度的变化在逐渐减小血液由弯管的内侧流向外侧，流动轨迹由曲线渐渐变为直线。图5-7为在雷诺数为760时U型血管150度时的压强横截面等高图，此时的压强梯度已不再剧烈变化，只是有些轻微的变动。我们再来看当血液流出U型血管出口处直管压强的等高图也就是图5-8，180度时的压强横截面等高图，我们会发现此时的压强梯度趋于稳定。 6.结论通过对七幅图的观察，我们不难发现当血液刚刚进入和流出U型血管时，压强梯度基本趋于恒定，变化很小，而当血液进入弯管后和流出弯管时，由于其运动轨迹发生变化由直线变为曲线和由曲线变为直线，其压强梯度变化尤为剧烈，此时最容易形成二次流。而由于二次流的作用，血液会和血管壁进行多次碰撞，容易发生回流，大分子颗粒容易发生堆积，从而造成血管壁增生，进而为医学提供参考。总结与体会时光的流逝也许是客观的，然而流逝的快慢却纯是一种主观的感受。当自己终于可以从考研、找工作、毕业论文的压力下解脱出来，长长地吁出一口气时，我忽然间才意识到，原来四年已经过去，到了该告别的时候了。一念至此，竟有些恍惚，所谓白驹过隙、百代过客云云，想来便是这般惆怅了。可是怅然之后，总要说些什么。大学四年，生活其实很简单，只是一些读书、写字和考试的周而复始。如果把这种单调的生活看作一场场循环的演出，提供那么我只是一个安静的演员。这篇毕业论文也称不上什么精彩的台词，只不过是这种循环演出即将告一段落时的谢幕词。但是无论多么蹩脚的演员，无论台下有多少观众，即使是只说给自己听，在他谢幕时也总要感激一些人，是这些人帮助他走上舞台，成功或者不那么成功地“演出”。谢 辞我在这里首先要感谢的是我的毕业论文指导老师老师。这篇毕业论文从开题、资料查找、修改到最后定稿，如果没有她的心血，尚不知以何等糟糕的面目出现。我很自豪有这样一位老师，她值得我感激和尊敬。 感谢和我共度四年美好大学生活的2009级应用物理专业的全体同学。感谢理化学院的所有授课老师，你们使我终身受益。感谢所有关心、鼓励、支持我的家人、亲戚和朋友。参考文献1 Chen S, Doolen G D. Lattice Boltzmann method for fluid flows J. Annual Review of Fluid Mechanics, 1998, 30(1): 329-364.2 Qian Y H, Deng Y F. A lattice BGK model for viscoelastic mediaJ. Physical review letters, 1997, 79(14): 2742-2745. 3 Ladd A, AJC. Numerical simulations of particulate suspensions via a discretized Boltzmann equationJ. Journal of Fluid Mechanics Digital Archive, 2006, 271: 285-309.4 Swift M R, Orlandini E, Osborn W R, et al. Lattice Boltzmann simulations of liquid-gas and binary fluid systemsJ. Physical Review E, 1996, 54(5): 5041-5052.5 Dawson S P, Chen S, Doolen G D. Lattice Boltzmann computations for reaction diffusion equationsJ. The Journal of Chemical Physics, 1993, 98: 1514 6 Rakotomalala N, Salin D, Watzky P. Simulations of viscous flows of complex fluids with a Bhatnagar, Gross, and Krook lattice gasJ. Physics of Fluids, 1996, 8: 3200.7 Gabbanelli S, Drazer G, Koplik J. Lattice Boltzmann method for non-Newtonian (power-law) fluidsJ. Prog. Physical review. E, 2003, 72: 46312.8 Watanabe T. Flow pattern and heat transfer rate in RayleighBnard convectionJ. Physics of Fluids, 2004, 16: 972.9 Chen H, Kandasamy S, Orszag S, et al. Extended Boltzmann kinetic equation for turbulent flowsZ. American Association for the Advancement of Science, 2003. 10 Mit