计算流体动力学仿真（精选8篇）
计算流体动力学仿真 篇1
摘要：高炉炼铁过程是我国冶金工业中的重要组成部分, 在炼铁过程中包括了多种物理变化与化学反应, 并且在高炉炼铁过程中多种物质不断相互作用, 所以高炉被化工领域称作是最复杂的反应器之一。当前高炉内部现象的控制研究工作对优化高炉炼铁工艺具有重要意义, 因此引起了众多专家的关注。然而随着高炉规模的逐渐扩大, 内部的动态难以掌握, 不能为专业人员提供相关信息。该文通过计算流体力学角度来解析高炉炼铁过程, 提出对我国现阶段高炉炼铁建模仿真框架, 对数据可视化处理, 还原高炉内部状况。
关键词：高炉炼铁,控制研究,计算流体力学,建模仿真
随着社会的发展, 钢铁行业的竞争也愈来愈明显, 许多欧洲国家开始自上而下地产业结构调整, 纷纷淘汰落后的工艺与系统, 试图最大范围地降低冶炼成本, 提升企业竞争力。从二十世纪七十年代以来, 许多欧洲国家纷纷建设起两千立方米以上的炼铁高炉, 我国也开始调整钢铁产业的发展政策, 认为应当淘汰1000立方米以下的高炉, 钢铁产业新建的高炉必须是大高炉。随着高炉内部的传感器的发展与应用, 让我们在对高炉正确使用过程中可以了解到内部的动态特性[1]。然而由于冶炼过程中, 高炉内部具有过高的温度, 因此传感器难以有效发挥出自身的效果, 因此需要我们进行建模对其状态进行研究
1 高炉过程的建模理论
1.1 炼铁过程概述
高炉炼铁过程包括了传热、传质以及化学反应等多种变化与过程, 涉及到热力学、化学、物理、传质等专业理论。高炉炼铁是将开采的铁矿石还原为生铁的过程, 人们将铁矿石、焦炭等资源按照科学配比后从高炉顶端的装置送到高炉内部, 高炉中的矿石与焦炭会出现分层的结构变化。然后铁矿石会在下降过程中逐渐产生化学反应被还原为铁与废弃物质, 最终聚集在高炉底部, 分别从铁口与渣口运送出来[2]。 (如图1)
从实质上来分析高炉冶炼过程就是在一定条件与环境下, 将铁从矿石原料中提炼出来, 这种冶炼技术是从古人发明的方法中总结与改进的, 具有生产量大与便捷的特点, 是当前钢铁行业最主要的冶炼方法, 每年生产量占全世界生产总数的95%以上。
1.2 矿石还原反应模型
我们冶炼的高炉从根本上来说是一种特殊的化学反应容器, 冶炼的核心就是让化学反应顺利发生。一般产生的化学反应分成三种:物质还原反应、C气化反应与水煤气变换反应, 对于模型仿真中的化学反应我们将计算过程主要集中在物质计算化学的反应速率方面。为了能够研究出矿石中铁元素的还原原理, 我们可以通过三界面未反应核模型理论进行冶炼过程[5]。 (如图2)
这个原理是在还原反应中进行的, 随着化学反应的逐渐深入, 铁元素没有反应的部分会朝着核中心缩小, 在最后形成了不会反应的核心。根据氧化铁逐层反应的特征, 我们可以看到一个矿石球在反应到一定程度后就会逐渐变为明显的三个界面, 分成四层, 分别为Fe-Fe O-Fe3O4-Fe2O3。而在这个核心最外层就是化学反应之后的最后产物, 具有多孔的Fe层物质。
我们可以将还原过程分成下面的环节:将还原的气体扩散、对内扩散、吸入、进行化学反应、氧化后气体吸收分解、对内扩散、外扩散等, 而整个过程中最慢的环节会决定所有化学反应中的速率, 也就是我们所说的速度控制步骤:矿石周围气体层的扩散运动、内部多孔产物的扩散作用与表层化学反应过程。
1.3 化学反应模型的实现
我们要想在冶炼中获得效果较好的化学反应, 就需要解决许多问题, 比如高炉中的化学反应模型并不是独立的, 而是会受到高炉内部许多因素的影响, 比如物质密度、温度、速率等方面, 而我们如果控制好, 就可以反过来影响这些因素。并且我们为了让高炉模型收敛, 就需要认识到不同反应是在哪种情况下放生的。
当内部温度低于848℃后, Fe3O4会被内部的氢气或者一氧化碳分解还原为铁;当内部高于848℃后, 就会经过变化为Fe O的过程。当Fe O在冶炼过程中的还原几率小于1的时候, 就会进行界面反应;反之, 就会进行二界面反应。
1.4 高炉过程的多相流关系
在近几年, 我国高炉炼铁在测量方式、技术等方面得到了巨大的概述, 所以需要构建多相态的数学模型才可以更加有效地了解到高炉内部的情况。由于内部流动现象复杂多变, 一般人们会分成软熔带、回旋区与死区来代表不同区域的变化, 让高炉的变化过程容易理解起来。然而在冶炼过程中的一些变化就会引起整个过程的结果变化, 所以需要进行整体的多相流分析。一般炼铁分成了四种相态, 固态、液态、粉态与气态。而根据四种相态要根据边界条件来合理约束闭合模型, 把握好各个相态之间的耦合关系, 进而顺利地实现的化学反应, 为高炉模型的仿真打下稳固的基础。
1.5 相态间的动量传递
高炉内部, 焦炭与矿石运送到高炉内不, 会形成交替的层结构, 而气流通过填料床内部, 会受到不同方向的曳力。高炉内部的矿石与焦炭形成的颗粒大小不同, 因此产生的曳力也不同。而在高炉下不, 残渣液体会随着温度升高而落下, 这样上升的煤气流会与残渣液体出现明显的耦合作用, 产生了液态与气态的动量交换。
2 高炉模型的仿真过程
当气体从鼓风口进来之后, 气体会受到高炉内软熔带的影响而移至高炉内部中轴线的方向, 导致中部出现了气流背离炉壁的现象。而受到炉内倾斜的状态与气体中较大的固体颗粒会导致孔隙度增大, 让中轴线地方的气流较强。 (如图3)
如图4所示, 炉内三个界面的各个反应率不同, 随着内部固体的下降, 逐渐会吸收到来自高炉气体的温度, 会与还原性的气体相接触, 接着产生了还原反应。
当固体从顶部进入到炉内, 随着高炉内部的半径逐渐增加, 因而固体的速度也会慢下来。在靠近高炉的出口处, 逐渐会增加速度。由于冶炼的高炉内部与死区表面会有相互作用力, 因此在出口处, 固体的方向会指向出口。图5表示为程序循环的100次、200次与500次, 即最终稳态固体的速度矢量的结构[4]。我们在实际测量速度的时候会采用TDMA的方式, 而数据从边界会传递到内部, 这导致顶端的速度变化程度较为剧烈。
由于固体在向下运动的时候会不断吸收气体中的热量, 所以会导致温度升高, 产生一系列化学物理变化。由于我们在参数设置方面会与实际设置有一些差别, 因此对于更精确的数值还需要给模型框架更准确的初值才可以计算出来。
3 结束语
该文通过计算流体力学的理论来分析高炉炼铁过程, 通过质量、动量、能量过程分析, 考虑到实际高炉内部产生的物理与化学变化, 对高炉炼铁过程建模仿真框架构建。这种方法能够成功地预测出高炉炼铁内部的化学反应状况与运动过程, 并研究出一些操作变量对整个高炉情况的影响, 能够客观科学地反映出实际工作情况, 对我国高炉自动化实现打下了良好的基础。
参考文献
[1]高可政, 马忠亮, 萧忠良.变燃速发射药芯料体积流率波动值的数值模拟[J].火炸药学报, 2010 (1) .
[2]薛文良, 魏孟媛, 陈革, 程隆棣.喷气织机主喷嘴内部流场的数值计算[J].纺织学报, 2010 (4) .
[3]邹惠芬, 张培红, 叶盛.流体力学多媒体教学的探讨[J].沈阳建筑大学学报:社会科学版, 2008 (4) .
[4]韩向科, 钱若军.流体力学基本理论研究综述[J].空间结构, 2008 (3) .
[5]杨树人, 王春生, 杨英.黏弹性流体在内管做轴向运动的偏心环空中的速度分布[J].大庆石油学院学报, 2005 (1) .
[6]肖进新, 田翼, 顾永涛, 罗妙宣.应用减压型毛细管粘度计研究表面活性剂溶液的流变性质[J].实验室研究与探索, 2000 (1) .
[7]韩式方.非牛顿钻井流体流动计算机智能解析方法[J].西南石油学院学报, 2001 (4) .
[8]白鹏, 钱庆利, 王红星, 袁绍军.水蒸气精馏流体力学和传质的实验研究[C]//第三届全国传质与分离工程学术会议论文集, 2002.
投身计算流体力学 篇2
在这里，他投身于计算流体力学的研究及教学中，并对湍流这一百年难题进行了深入探索。
致力于计算流体力学软件开发
计算流体力学（Computational Fluid Dynamics，CFD），是采用数值计算的方法研究流动问题的一门学科。该学科在工程及科研中发挥了重要作用，尤其是航空航天领域。
飞机、导弹等飞行器设计中需要解决大量的空气动力学问题，早期主要依赖风洞试验，随着计算流体力学及超级计算机的产生及发展，采用数值计算的方法解决这些气动问题成为可能，这就是计算流体力学（CFD）。如今，CFD已在飞行器设计中发挥巨大的作用，很多飞行器设计仅需要少量风洞试验即可完成定型。CFD计算需对软件需求强劲，目前国际上已经有了很多较为成熟的商业CFD软件，而像波音公司这样的大型企业也有自己的内部CFD软件。与国外比，我国在CFD软件开发及推广方面仍有一定差距。
面向应用中的强劲需求，李新亮在CFD算法、模型及软件方面投入十多年的研究，开发了一套开源的计算流体力学软件Open CFD。该软件包含了有限差分及有限体积两个求解器。
有限差分求解器的最大特点是具有非常高的精度，最高精度可达十阶。该软的差分库包含了其自行开发的多种高精度差分格式，以及目前流行的高精度差分格式。主要用于直接数值模拟及大涡模拟等复杂流动的机理研究。有限体积求解器采用多块结构网格，集成了目前常用的数值方法及湍流模型，可用于像整架飞机这样的复杂外形流动的数值模拟。
该软件的另外一个特点是，具有很强的并行可扩展性。测试结果显示，该软件在常规CPU体系上实现了五万余CPU核心的高性能并行计算；而使用众核系统测试时，并行规模达到近百万CPU核心。该软件获得了中科院网络中心2011年度的“计算规模奖”和2013年度“最佳应用奖”。
该软件向国内外开放源代码，二维代码可在流体中文网下载（http：//www.cfluid.com/bbs/forumdisplay.php？fid=2），目前下载量达到上万余次。三维代码也被国内外20多家单位使用。
专注CFD研究和教学
李新亮及其团队在CFD算法，特别是在高精度激波捕捉格式方面取得了丰富的研究成果，先后开发了优化保单调差分格式、加权群速度控制格式等一系列高精度激波捕捉方法。与国际上流行的高精度方法相比，这些方法在尺度分辨率及计算效率方面具有突出表现。
李新亮的精力不仅花在科研上，在教学方面，他更是以自己的一言一行深深影响周围的每一位学生，他教导学生就像一个农民侍弄自己的庄稼一样，用爱心，用耐心，教他们为学、为人。
爱因斯坦曾说，对一切来说，只有兴趣才是最好的老师。兴趣的培养和激发是老师的责任和义务。在这方面，李新亮很用心。每年春季，他在力学所开设计算流体力学的课程，总是吸引了周围高校和研究所的很多学生前来旁听，外校听课学生的比例甚至超过了本单位的学生。他的课件和录像在“流体中文网”等网站被下载上万次，成为学习《计算流体力学》的主流课件。
探索湍流百年难题
湍流是流体的一种多尺度复杂运动状态，自然界和工程中的流动多为湍流状态。湍流是自然科学中的百年难题。到目前为止，湍流的形成、发展机理及其数学模型仍未完全明确，该问题的研究一直是流通力学的学科前沿。飞机、导弹等飞行器设计过程中，无论是气动外形还是发动机设计都会遇到大量与湍流（包括转捩）有关的问题。这些湍流问题是飞行器设计中的“卡脖子”问题，也是工程关注的重点问题。
面向湍流难题，李新亮利用自主开发软件结合大规模并行计算，采用直接数值模拟（DNS）手段进行了深入研究。直接数值模拟是采用密集的计算网格计算出湍流全部尺度流动细节。该方法不依赖湍流模型，因而可以准确地计算湍流，且可以给出湍流的全部信息。是研究湍流机理、实现控制及及建立和改进湍流模型的有效手段。由于湍流的多尺度性，分辨湍流的全部细节需要非常密度的计算网格（网格点数动辄数亿，甚至更多）因而计算量非常巨大，对计算资源及计算软件要求苛刻。而飞行器湍流属于可压缩湍流、包含了湍流与激波等多种复杂因素，其直接数值模拟难度更大。国际上该流动的DNS研究也是近20年内的事情。
面向飞行器可压缩湍流的难题，李新亮及其团队先后对可压缩槽道、平板、钝锥、压缩折角、激波—湍流边界层干扰等复杂问题进行了直接数值模拟研究。在此基础上建立了可压缩湍流数据库，该湍流数据库为湍流机理及湍流模型研究奠定了条件。到目前为止，该湍流数据库为国内外多个科研小组提供了研究可压缩湍流的第一手数据。
钝锥体是火箭及导弹等高速飞行器的典型头部外形，其转捩位置及转捩机理在航天领域倍受关注。直接数值模拟（DNS）是研究该问题的有效手段，但由于计算复杂，国际上尚没有高超声速有攻角钝锥边界层湍流DNS结果报道。面向国家需求，李新亮团队开展了来流马赫数6小攻角球钝锥边界层转捩到湍流的直接数值模拟，得到了整个钝锥边界层转捩过程的全部流场。发现了小攻角情况下，高速钝锥表面上的转捩位置呈现出非单调曲线分布。在此之前，尚没有高超音速钝锥边界层湍流直接数值模拟（DNS）的公开报道，该结果是首个公开报道的高速钝锥边界层湍流DNS算例。该研究为认识小攻角高速钝锥转的捩机理提供了有力依据。
在求索中创新，收获中进步
科学研究是探求未知的工作，贵在求索创新。在投注大量科研精力的领域中，李新亮一路耕耘一路收获。他承担了自然科学基金、“973”“863”等课题以及国防专项等一系列科研项目。获得2000年度中国科学院自然科学奖二等奖、2013年陕西省国防科技进步一等奖等奖项。
除了收获这些荣誉，他还在多项社会活动中发挥所长。现在，李新亮兼任中国空气动力学会物理力学专业委员会副主任委员，中科院超级计算用户委员会委员，“十二五”中科院信息化专题用户组成员，计算物理学会常务理事，《计算物理》杂志编委等众多职务。目前还担任国际知名杂志“Computers & Fluids”的地区编委。
计算流体动力学仿真 篇3
1 建模
非接触式承片台一般都采用伯努利原理来实现,对于理想不可压缩的流体,存在以下方程:
式中:ρ-流体的密度;υ-流体的速度;g-重力加速度;h-流体处于的高度;p-流体所受的压强。
常温常压下的空气可视为理想不可压缩气体,并忽略空气的重力,则上式可简化为
上式左侧两项分别称为动压和静压,即沿流线运动的理想气体的动压和静压的总和为常数。
本文所述的非接触式承片台可用于放置12英寸的晶圆,其外径尺寸为340mm,其结构如图1所示。
该承片台具有一个进气口,直径为8mm,布置在承片台下部,采用中空轴电机供气。承片台腔体内部共有6个出气口,呈圆周均匀排布,直径为3mm,通过内部管路将气体由底部送入腔体内部。顶端面上圆周均布4个圆柱销,防止晶圆的径向移动。当晶圆放置在承片台上时,气流在承片台和晶圆所形成的空腔内部高速旋转,腔内空气动压增大,静压减少,最终气体通过晶圆与承片台之间的间隙排出。因为晶圆下表面的静压小于上表面的静压,这样就实现了将晶圆非接触式地压在承片台上方的目的。
2 划分网格
模型网格的划分是流体动力学分析中十分重要的一步,因为网格质量直接关系到计算结果的准确性和收敛性。网格的划分可分为结构网格和非结构网格两种。结构化网格是指网格区域内所有的内部节点都具有相同的毗邻单元,即与网格区域内的不同节点相连的单元数目都相同[2]。它的主要优点有:(1)可以很容易地实现区域的边界拟合,适于流体和表面应力集中等方面的计算;(2)网格生成的速度快;(3)网格生成的质量好,数量少;(4)数据结构简单,计算收敛速度快;(5)对曲面或空间的拟合大多数采用参数化或样条插值的方法得到,区域光滑,与实际的模型更容易接近。但它相应的缺点就是对复杂模型的划分比较困难。所以很多人为了方便直接采用四面体单元自动划分,但在处理壁面附近层时,四面体单元会在某一方向上被拉长,导致网格质量变差,方程不易收敛,为了提高网格质量,只能使用更小尺寸的四面体网格,但这将会带来单元数目呈几何倍数的增长,增加了计算时间。因此建立六面体网格会明显地提高计算的收敛速度和准确性。在此对模型进行结构化网格划分。
为了进一步减少网格数量,节省计算时间,这里使用ANSYS ICEM模块对模型的一半进行网格划分。步骤如下:首先建立一个Block并将它映射到上部的半圆柱上,然后进行O-grid的划分。使用Extrude Face命令将划分好O-grid的中心面向下拉伸出一个新的Block,并将新的Block映射到下部的圆柱面上。对新Block进行O-grid划分,这就完成了主体部分的结构网格划分,因为入口处管路直径很小只有3mm,所以此处网格需要更细致的划分。用Split Block命令将入口管路附近处的Block切割出一个和管口大小相近的小Block,再次使用Extrude Face命令将小Block的面向外拉伸出一新Block。将此Block映射到入口管路上,调整Vertex点位置,再进行O-grid划分,即完成了入口处的网格划分,如图2所示。用同样的方法完成其他两处入口网格,并对壁面处网格进行了加密处理。最终完成网格划分,如图3所示。
3 仿真结果分析
对划分网格的模型设定相关边界条件以后就可以进行仿真计算了。网格模型的入口是物理模型的出口,即6个直径为3mm的管路,相对压力0.2MPa,晶圆与承片台之间的间隙即为出口,相对压力0MPa。其它壁面均设置成No slip wall。最终计算晶圆下表面压力云图如图4所示。
由图4可以看出,晶圆的边缘处,也就是气体的出口部分相对压强为正,而在腔体内部的部分相对压强为负,整个表面压强呈梯度分布,越靠近中心负压值越大,最大相对负压约为-86Pa,将负压值乘以相应区域的面积并求和,即可得到晶圆上表面所受到的压力值。
通过调整出口间隙值、腔体直径、腔体深度等参数,重新建立模型并求解,最终可得到一组较优的参数值,使在进气流量相同的情况下,晶圆获得最大的下压力。此时出口间隙为2mm,腔体直径为280mm,腔体深度为22mm。
4 实验验证
为了验证仿真结果的准确性,将对实际生产出的承片台进行试验验证,实验环境如图5所示。
实验仪器:数显吊钩秤,量程2kg,精度1g。
实验方法:(1)用细线
将晶圆吊在数显吊钩秤的正下方,注意保持晶圆的水平,不要偏斜;(2)将晶圆放置在承片台上并打开进气开关,仔细调整压力阀使进气压力达到0.2MPa,此时晶圆被吸附在承片台上方;(3)手提着数显吊钩秤缓慢提起,直到晶圆在向上拉力作用下完全脱离承片台,记录显示屏上出现过的最大拉力值。向上提起过程中要注意保持数显吊钩秤的垂直。
因为检测方法比较简单,为减少人为误差共检测3次并取其平均值,数值见表1。
由竖直方向力的平衡可知,数显吊钩秤的最大拉力等于晶圆的重力和承片台的吸附力之和。已知晶圆自重为126g,可得吸附力为405.3-126=279.3g=2.74N。由仿真结果计算出的吸附力约为3.44N,仿真值比实际值大了25.5%。分析误差主要是由人为的读数误差、晶圆抬起过程中的偏心受力和进气口的压力误差等引起的。实际吸附力值达到了晶圆自重的2倍以上,可以满足使用要求。
5 结论
(1)通过仿真结果的压强云图可以发现,承片台中心区域为负压,最小值为-86Pa,可以使晶圆有效地被夹持,同时晶圆的边缘为正压,可以有效地阻止外部杂质进入,从而保证了晶圆的背面不被污染。仿真结果与实验结果相比较为接近,与传统设计方法相比提高了设计效率,降低了生产成本,最终产品达到了预期的效果。
(2)ANSYS中ICEM模块强大的六面体网格编辑功能,可使模型网格达到更高的质量,让仿真结果更加接近真实值,提高准确度。
摘要：介绍了一种伯努利原理非接触式承片台设计的新方法,即使用流体动力学软件对承片台的结构和性能进行分析,并详细叙述了承片台结构网格的划分方法。根据仿真结果对其结构进行改进设计,设计产品最终达到了良好的使用效果。
关键词：CFD技术,伯努利承片台,设计方法
参考文献
[1]王福军.计算流体力学分析[M].北京:清华大学出版社,2004.
[2]万映娟.离心泵叶轮内流场的数值模拟[D].成都:西华大学,2008.
[3]徐蕾,罗会信.基于ANSYS ICEM CFD和CFX数值仿真技术[J].机械工程师,2008(12):65-66.
计算流体力学发展综述 篇4
计算流体力学 (Computational Fluid Dynamics, 简称CFD) 是当代迅速发展的一门学科, 是利用高速计算机求解流体流动的偏微分方程组, 目的是为了更好的从定性上和定量上了解流体流动的物理现象, 改进设计的一门学科。目前在航空航天、交通运输、造船、气象、海洋、水利、液压和石油化工等工程领域都有广泛的应用。
作为流体力学的新手段, CFD是一种令人鼓舞的模拟流体流动的方法, 它大大缩短了设计的时间, 节省了设计费用。它相对于理论方法来说, 具有假设限制少、应用范围广的特点, 其方法也容易应用。相对于实验来说, 计算流体力学很少有马赫数和物体尺寸的限制, 并且具有较高的经济价值。数值仿真优于实验的地方还在于:计算机仿真的诊断“探测”并不干扰流动且不使所研究的现象变得不可捉摸。CFD已经代替了许多环境列车的试验项目, 而试验的目的也逐渐从验证设计参数的合理性, 改变为对CFD数值模拟的正确性及最终设计的校核。CFD不仅可以为列车提供快速而经济的设计依据, 并且可以观测到一些试验中无法观测到的物理现象, 还可以为高速列车的设计提供理论依据。
计算流体力学是在飞机工业的需要的基础上发展起来的。由于要求用流体力学理论来了解和指导飞机设计, 当时, 由于飞行速度很低, 可以忽略粘性和旋涡, 因此流动的模型为Laplace方程, 研究工作的重点是椭圆型方程的数值解。利用复变函数理论和解的迭加方法来求解析解。随着飞机外形设计越来越复杂, 出现了求解奇异边界积分方程的方法。以后, 为了考虑粘性效应, 有了边界层方程的数值计算方法, 并发展成以位势方程为外流方程, 与内流边界层方程相结合, 通过迭代求解粘性干扰流场的计算方法。
同一时期, 许多数学家研究了偏微分方程的数学理论, Hadamard, Courant, Fried richs等人研究了偏微分方程的基本特性、数学提法的适定性、物理波的传播特性等问题, 发展了双曲型偏微分方程理论。以后, Cou-rant, Friedrichs, Lewy等人发表了经典论文, 证明了连续的椭圆型、抛物型和双曲型方程组解的存在性和唯一性定理, 他们还给出了著名的稳定性判别条件:CFL条件。Von Neumann, Richtmyer, Hopf, Lax和其他一些学者建立了非线性双曲型方程守恒定律的数值方法理论, 为含有激波的气体流动数值模拟打下了理论基础。
随着技术的发展, 仅采用当时流体力学的方法, 研究较复杂的非线性流动现象是不够的, 特别是不能满足高速发展起来的宇航飞行器绕流流场特性研究的需要。针对这种情况, 一些学者开始将基于双曲型方程数学理论基础的时问相关方法用于求解宇航飞行器的气体定常绕流流场问题, 这种方法虽然要求花费更多的计算机时, 但因数学提法适定, 又有较好的理论基础, 且能模拟流体运动的非定常过程, 所以在60年代这是应用范围较广的一般方法。以后由Lax、Kreiss和其他著者给出的非定常偏微分方程差分逼近的稳定性理论, 进一步促进了时间相关方法。当时还出现了一些针对具体问题发展起来的特殊算法。
我国在20世纪50年代也开始了计算流体力学方面的研究。我国早期的工作是研究钝头体超声速无粘绕流流场的数值解方法, 研究钝头体绕流数值解的反方法和正方法。以后, 随着我国宇航事业的发展, 超声速、高超声速绕流数值计算方法的研究工作发展很快。对定常欧拉方程数值解的计算方法进行研究, 并给出了钝体超声速三维无粘绕流流场的计算结果。
20世纪70年代, 在计算流体力学中取得较大成功的是飞行器跨音速绕流数值计算方法的研究。首先是Murman和Cole用松弛方法求解位势流小扰动方程, 数值模拟带激波的跨声速绕流场。解决了跨声速绕流中的混合问题。在他们的工作中第一次将迎风格式应用于空气动力学问题的模拟。不久以后Jameson提出了旋转格式, 将穆尔曼一科勒方法推广于求解三维跨声速绕流的全位势流方程, 获得成功。同一时期, 我国开展了采用时间相关方法求解非定常欧拉方程、可压缩N-S方程和简化N-S方程的计算方法研究。
在差分格式的构造方面, 提出了求解欧拉方程的特征符号分裂法和三层格式等。在可压缩N-S方程的求解中, 计算方法有了很大进展, 先后提出了开关函数法、调解因子方法、紧致迎风格式、推进迭代法、无波动无自由参数的耗散格式、界值为限格式和耗散比拟方法等。这些研究工作进一步改进了计算方法精度, 提高了求解效率, 且对流场激波的数值模拟有较高的分辨能力。而且这些研究成果使得我们在计算流体力学的差分方法研究工作中初步形成了自己的特点。
2 计算流体力学的现状
进人20世纪80年代以后, 计算机硬件技术有了突飞猛进的发展, 千万次机、亿次机逐渐进人人们的实践活动范围。随着计算方法的不断改进和数值分析理论的发展高精度势值模拟已不再是天方夜谭。同时随着人类生产实践活动的不断发展, 科学技术的日新月异, 一大批高新技术产业对计算流体力学提出了新的要求, 同时也为计算流体力学的发展提供了新的机遇。实践与理论的不断互动, 形成计算流体力学的新热点、新动力, 从而推动计算流体力学不断向前发展。
首先, 在计算模型方面, 又提出了一些新的模型, 如新的大涡模拟模型、考虑壁面曲率等效应的新的湍流模式、新的多相流模式、新的飞行器气动分析与热结构的一体化模型等。这就使得计算流体力学的计算模型由最初的Euler和N-S方程, 扩展到包括湍流、两相流、化学非平衡、太阳风等问题研究模型在内的多个模型。其中以考虑更多流动机制, 如各向异性的非线性 (应力/应变关系) 湍流研究为重点。研究结果再次证明, 万能的湍流模型还不存在, 重要的是如何在模型精度和计算量上较好地取得折中;也有学者从更高层次研究湍流模型问题, 由湍流流动中速度不可微, 怀疑N-S方程的有效性, 进而提出以积分方程为基础的数学模型。
其次, 在计算方法方面, 又提出了一些新的计算方法, 如新的遗传算法、无网格算法、新型高精度紧致格式、气动计算的新变分原理、结构/非结构混合网格新技术、新型动网格技术等等。目前计算方法的研究集中在高精度格式方法, 即追求三阶精度以上, 其中又以解决真正实际问题。除此之外, 计算方法研究还涉及带限制器的高阶插值、谱方法、拉格朗日方法, 时空守恒元方法等等。将其它方法引进传统的计算流体力学也是现阶段的重要成果之一, 其中特别值得一提的是将基因算法与传统计算流体力学结合在一起, 在域分裂和最优化设计等许多方面显示出了良好的应用前景。在算法分析上, 除传统的精度、稳定性、收敛性等方面的分析, 还有更深层次的数值动力学分析, 即将数值方法看成是动力系统来进行分析, 揭示了许多奇异的数值现象。
再次, 在研究成果方面, 英国M.A.Lesdhziaer关于湍流模型、美国H.C.Yee关于计算不确定性、日本学者的玻耳兹曼方程解流动问题、德国的E.von Lavante关于使用并行计算机进行发动机气缸流场涡和激波的非定常流动模拟等等, 都有较新的学术思想, 较高的学术水平。
目前, 计算流体力学研究的热点是:研究计算方法, 包括并行算法和各种新型算法;研究涡运动和湍流, 包括可压和不可压湍流的直接数值模拟、大涡模拟和湍流机理;研究网格生成技术及计算机优化设计;研究计算流体力学用于解决实际流动问题, 包括计算生物力学、计算声学、微型机械流动、多相流及涡轮机械流动的数值模拟等。
3 CFD软件技术
CFD技术艰深的理论背景与流体力学问题的复杂多变阻碍了它向工业界推广。一般工程技术人员很难较深入地了解这门学科, 由专家编制的程序用起来也不容易, 因为总有不同条件、参数要根据具体问题以及运算过程随时做出修改调整, 若不熟悉广法和程序, 往往会束手无策, 此外, 前、后处理也显得十分棘手。CFD研究成果与实际应用的结合成为极大难题, 这一切曾使人们对CFD的工程应用前景产生疑虑。在此情况下, 通用软件包应运而生, 使CFD计算变得方便、简单。
CFD软件一般包括三个主要部分:前处理器 (建模, 网格生成等) , 解算器 (具体的数值运算) 和后处理器 (运算结果的具体演示) 。常见的CFD软件有:FLUENT, PHOENICS, CFX, STAR-CD, FIDAP等。
以FLUENT公司开发的大型CFD软件FLUENT为例, 它可计算从不可压缩 (低亚音速) 到轻度可压缩 (跨音速) 直达高度可压缩 (超音速) 流体的复杂流动问题。FLUENT本身所带的物理模型可以准确地预测层流、过渡流和湍流多种方式的传热和传质, 化学反应, 多相流和其它复杂现象。它可以灵活地产生非结构网格, 以适应复杂结构, 并且能根据初步计算结果调速网格。前处理软件Gambit提供了多方位的几何输入接口。计算采用有限容积法。通过图形后处理软件, 可以得到二维和三维图象, 包括速度矢量图、等值线图 (流线图、等压线图) 、等值面图 (等温面和等马赫面图) 、流动轨迹图, 并具有积分功能, 可以求得力和流量等。
4 展望及结论
计算流体力学主要向两个方面发展:一方面是研究流动非定常稳定特性、分叉解及湍流流动的机理, 更为复杂的非定常、多尺度的流动特征, 高精度、高分辨率的计算方法和并行算法;另一方面是将计算流体力学直接用于模拟各种实际流动, 解决工业生产中提出来的各种问题。美国和日本在这两方面做得最为突出。在我国经济飞速发展的今天, 一些计算流体力学问题的解决, 将有利于我国的国民经济建设工作, 我们需要迎头追赶。
参考文献
[1]傅德薰.第二届国际计算流体力学会议简介[J].力学进展2003 (3) :388.
[2]刘国俊.计算流体力学的地位、发展情况和发展趋势[J].航空计算技术, 1994 (1) :15-21.
[3]魏淑贤等.计算流体力学的发展及应用河北理工大学2005, 27 (2) :115-122.
[4]王东屏等.CFD数值仿真在高速列车中的应用及验证2005Fluent中国用户大会论文集.
计算流体动力学仿真 篇5
《计算流体动力学及其应用》是能源动力方向硕士和博士研究生的一门学科基础课程。自上世纪90年代开设以来, 前后经历了三次大的变革。在1995—2000年期间, 课程名为《叶轮机械的数值模拟计算》, 只针对能源动力学院流体机械专业方向的研究生开设;自2001年起, 课程改名为《计算流体动力学及其应用》, 面向能源动力一级学科及其下属各二级学科的研究生, 并成为能源动力学科方向研究生的一门基础课程;2006年以前, 课程授课内容以计算流体动力学方面的原理和方法为主, 选课学生主要为能源动力学科方向的研究生;从2006年开始, 为适应广大研究生的选课需要, 我们对教学内容进行了适当调整, 辅以CFD商用软件的实例和应用, 以期实现理论性与实践性并重的教学理念, 并将课程面向全校研究生开放。随着计算机技术的飞速发展, 计算流体动力学的应用日益广泛。众所周知, 计算机硬件水平的提升, 将相应地促进CFD商用软件功能更加强大, 应用更加广泛, 最终使得CFD商用软件得到了前所未有的发展。同时, 随着研究生招生规模的扩大, 使得选修《计算流体动力学及其应用》课程的研究生人数大增, 从上个世纪90年代的十几个学生, 到现在的一百多个学生, 而且涉及众多学科, 比如船海、化工、建筑、电气、交通、材料、光电等。《计算流体动力学及其应用》课程的历史与现状在一定程度上给我们将要进行的教学改革提出了新的要求, 同时也为我们指明了课程建设的新方向, 值得我们深入思考, 并付之于实践。
二、课程定位
《计算流体动力学及其应用》作为一门研究生的学科基础课程, 我们在进行改革之前, 应该首先考虑它的定位。华中科技大学作为一所教育部的“985”和“211”的高校, 一直以“研究型”大学著称。学校对于研究生的培养非常重视, 导师为每一位研究生制定了详细的培养计划, 课程的选修均有所考量。基于选修《计算流体动力学及其应用》课程的研究生人数众多, 涉及的院系广泛, 经任课教师讨论, 申请学校研究生院同意, 决定将该课程定位为高水平研究生课程。所谓高水平研究生课程, 初步确立的含义为, 高水平的学者, 采用高水平的教材, 以先进灵活的形式教授课程, 旨在培养学生坚实宽广的理论基础和系统深入的专门知识。高水平课程在内容上应该具有基础性、专业性和前沿性, 前沿性可以体现在任课教师结合自己的科研实践, 在讲授中融入一些与课程相关的前沿内容。
三、教改实践
基于高水平研究生课程这样一个定位, 我们开始着手进行课程的教学改革, 具体内容包括:组建教学团队、改革教学内容、建设实践教学平台。首先, 组建一支高水平、高素质的教学团队。教学团队由三位教师组成, 他们均具有博士学位, 高级职称。其中, 课程负责人张师帅副教授, 长期从事计算流体动力学及其应用方面的教学及研究工作, 自2006年起, 一直担任该课程的主讲教师;任课教师郭照立教授, 是目前国内计算流体动力学方面的顶尖学者, 国家杰出青年基金获得者, 并具有较高的国际知名度。郭教授团队在国内外权威学术刊物和会议上发表科学论文100余篇, SCI收录90余篇, SCI引用1200余次;任课教师陈胜副教授是一位青年学者, 在格子Boltzmann算法研究方面颇有建树。将《计算流体动力学及其应用》课程建设成一门高水平研究生课程, 得到了教学团队中每一位教师的支持, 大家一致赞同经常开展教学交流, 学习先进的教学方法和教学手段, 进一步提高教学效果。其次, 我们对教学内容进行了改革。教学团队根据选课研究生人数众多, 涉及的学科方向广泛, 重新制定了《计算流体动力学及其应用》课程的教学大纲, 确保讲授内容的基础性、专业性和前沿性。课程主要内容包括:控制方程的离散化方法、流场的求解计算方法、湍流模型及其应用、网格生成与计算技术、复杂流动的介观模型和数值方法、格子Boltzmann算法及其应用、经典CFD软件的基本用法等。而对于控制方程的离散化方法, 将重点介绍有限差分法和有限体积法;对于流场的求解计算方法, 将重点介绍SIMPLE及其系列算法;对于湍流模型及其应用, 将重点介绍k-ε模型及其应用;对于网格生成与计算技术, 将重点介绍结构网格和非结构网格的生成方法以及并行计算方法。同时, 还将邀请国内外计算流体动力学方面的专家学者前来开展专题讲座。对《计算流体动力学及其应用》课程进行教学改革是全体任课教师的共同愿望, 大家积极讨论, 并提出在现有教材的基础上, 编写具有自己特色的教材等建议。在改革教学内容的同时, 教学团队还利用学院现有的条件, 建立“计算流体动力学”软件平台, 该平台拥有高性能的计算工作站, 可以开展并行计算、直接数值模拟等大型计算研究, 为研究生开展离散方法、网格生成方法、计算方法以及复杂流动模型等研究工作创造了良好条件, 同时也为对计算流体动力学方面的前沿研究课题感兴趣的大学本科生开展创新性研究工作提供了良好条件。与此同时, 该平台还拥有多种商用CFD软件, 比如FLUENT、CFX、STAR-CD、PHOENICS、Flo-EFD等, 成为广大研究生开展自主学习、自主实践、相互交流的优良场所。还可以根据研究生的需求, 安排任课教师不定期地通过软件平台为学生解惑答疑, 引导研究生探索创新, 提高学术水平。
众所周知, 研究生学术水平的高低是一所大学学术水平的反映, 更是一个国家科技创新能力的反映。研究生不仅需要扎实掌握专业基础知识, 更需要具有较强的创新意识和创新能力。目前, 高等学校在培养研究生创新能力、提高研究生学术水平方面还有待加强。为此, 本文提出了通过对《计算流体动力学及其应用》课程进行教学改革, 并将之建设成具有基础性、专业性和前沿性的高水平课程, 进一步培养研究生的创新能力, 提高研究生的学术水平。同时, 本文对实践过程中的一些具体措施和经验进行了探讨。
参考文献
[1]吴新颖.论研究生教育的课程设置[J].科技导报, 2004, (2) :40-42.
[2]李学艺, 钟佩思.机械类全日制专业学位研究生教育探索[J].中国电力教育, 2013, (1) :22-23.
[3]李静波, 柴育玲.研究生教育创新模式研究[J].教育理论与实践, 2010, (5) :10-11.
[4]别敦荣, 万卫.学术性学位与专业学位研究生培养模式的特性比较[J].研究生教育研究, 2011, (4) :77-80.
计算流体动力学仿真 篇6
近年来越来越多的研究人员开始投身于地下空间楼梯流的研究,已经取得了一些进展。早期,城市地下空间水灾的研究工作主要集中于实验研究,包括全尺寸阶梯流实验和缩尺寸实验[4,5,6,7],而数值模拟研究较少。真实实验研究对研究地下空间阶梯流动力学特征比较重要,但其实验设计及布置、实验数据收集及处理不仅费时而且比较昂贵,并且大尺度的实验还存在水灾隐患。近几年来,随着计算机技术和数值计算方法的快速发展,地下空间楼梯流的数值研究成为一个热点。在地下空间楼梯流数值模拟分析方面,主要借鉴阶梯溢流坝水流数值计算模型,因为地下空间水灾的水力特性与阶梯溢流坝水流比较类似。近年来,国内外基于网格的数值方法在阶梯溢流坝水流数值模拟上应用较广,如VOF方法耦合标准k-ε紊流模型、VOF方法耦合Realizable k-ε紊流模型、浅水波方程模型等[8,9,10,11]。对于适用于行人行走的中间坡度(阶梯坡度多为26°~30°)的楼梯水流的数值模拟主要借助上述模型,在地下空间洪水模拟中多数都对楼梯进行了简化。Yoneyama等[12]基于VOF方法对地下空间直行楼梯上洪水水流进行了模拟,计算得到的速度场比实验测量值偏小。莫伟丽[13]基于标准k-ε紊流模式并利用VOF方法捕捉自由表面构建了地铁站三维数值计算模型,在模拟楼梯与通道的直角转弯处变化剧烈的水流时存在较大误差。申若竹等[14]采用VOF方法构建了三维计算模型,对地下空间直行无休息平台的直行楼梯上洪水水流进行了模拟,数模计算的最大流速高于实验值。Shao等[15]利用FLU-ENT软件中VOF模型和Realizable k-ε紊流模型对不同倾角与不同复杂形态(直行、L型、有无休息平台等)楼梯上洪水水流特征进行了模拟分析,并探讨了有无休息平台的情况下楼梯上洪水水流对人员疏散的影响作用。
传统的基于网格或依赖网格的VOF等数值方法对楼梯流复杂的自由表面形态,如破波、水流与固壁碰撞后的飞溅融合等现象较难处理和模拟分析。鉴于此,近几年来有些学者尝试将比较适合处理大变形问题及存在移动交界面的问题无网格方法应用于模拟分析阶梯流。Gotoh等[16,17]首次将粒子法用于地下空间阶梯流的研究中,基于MPS方法构建了地下空间楼梯水流的数值模型,计算获得了三种典型的楼梯流流态,然而MPS方法半隐式直接对动量方程求解压力,每个时间步需要求解一个大的泊松方程组来获得各个粒子的压力值,计算很耗时。吴建松等[18]采用串行计算的SPH方法对小尺度阶梯流进行了模拟计算,获得了较好的阶梯流态,但该模型计算效率较低,且仍局限于计算尺度较小的阶梯流场景。SPH方法和MPS方法对计算区域与控制方程的离散近似原理相同,但求解压力方程时区别较大,SPH方法引入弱可压性概念通过显式求解状态方程获得各个粒子的压力值,比MPS方法计算效率要高。
为了提高地下空间楼梯流数值计算效率,本文采用SPH方法并借助GPU并行计算技术对地下空间楼梯洪水水流进行数值建模和计算,并着重对比分析串行SPH模型、并行SPH模型还有基于MPS粒子法的地下空间梯流模型的计算效率。
1 SPH基本理论
1.1 SPH基本原理
SPH方法的理论基础为积分插值,其基本方程的构造分为两个关键步骤:场函数核近似和粒子近似[19]。在SPH方法中,函数f(x)积分表达式由式(1)定义:
式(1)中,x为粒子位移,W(x-x',h)为光滑函数,h为决定光滑函数影响区域的光滑长度。
SPH方法研究的系统是由有限个具有一定体积和质量的粒子来离散进行表示的,由此引入了粒子近似,从而将由SPH核近似得到的连续积分表达式转化为支持域内所有粒子叠加求和的离散化形式。因此,式(1)经粒子近似可得到粒子i处的函数近似表达式:
同理,还可得到场函数导数的近似式:
式(3)中,ρj和mj分别为粒子的密度和质量;
1.2 SPH控制方程
流体力学中连续场的动量方程可表示为:
式(4)中,v为速度,g为重力加速度,p为压力,φ为扩散项。根据流体的黏性类型,扩散项可用不同的形式描述,相应的也有不同的动量方程描述形式,SPH方法中现在主要有三种形式:人工黏度模型(无黏流动)、层流黏性模型和湍流黏性模型。
由式(1)~式(3)中场函数核近似和粒子近似方法,可得层流黏性形式的SPH控制方程组:
连续性方程为:
动量方程为:
1.3 状态方程
标准SPH方法主要应用于处理可压缩流动问题。理论不可压流体实际上都是可压的,Monaghan[20]引入人工压缩率概念,构造了求解控制方程中压力的状态方程:
式(7)中:ρ0为初始参考密度,水流取ρ0=1 000 kg/m3,γ取介于1~7之间的常数,对水通常取值为7,B用于限制密度的最大改变量,由具体问题而定。
1.4 粒子位移的更新
粒子位移更新方程采用Monaghan[20]提出的XSPH方法:
式(8)中,ε为常数,介于0~1.0之间,Liu等[19]建议ε取0.3,该值适用于大多数情况。XSPH方法能使粒子的运动速度与相邻粒子的平均速度接近,能有效防止相互逼近粒子的穿透。
2 GPUSPH模型
SPH方法对处理处于具有复杂界面变形的水流具有很大的优势,但作为一种无网格粒子法,计算量比较大,每个时间步都要计算部分粒子的相互作用,这非常耗时。为了解决这一局限性,Hérault等[21]率先把基于CUDA的GPU高性能计算方法应用于SPH计算体系,构建了基于GPU并行计算的SPH计算模型。GPUSPH模型采用Nvidia公司开发的CUDA编程语言及C++和Open GL进行编程。具体的程序架构和实现细节请详见参考文献[16]。本文工作在GPUSPH模型基础上进行了改进,加入了楼梯建模的基元图形和改进了边界处理算法。
2.1 光滑核函数
GPUSPH模型中,光滑核函有三种选择:二次函数、三次样条函数、Wendland五次函数。通常情况下,SPH方法中积分近似的精度随着光滑函数的次幂增高越来越好,但其计算效率则反之。GomezGesteria等[22]对SPH方法中比较常用的几个光滑函数的特性进行了对比分析,最后指出Wendland五次函数在计算精度和计算效率上有一个很好的折衷,其表达式如下:
αd在二维和三维场景中的取值分别为7/(4πh2)和7/(8πh3)。
2.2 密度校正
GPUSPH模型中,采用了两种方法对密度进行校正:零阶Shepard校正和一阶MLS校正。本文的洪水模拟使用了一阶MLS校正,密度重新初始化过程采用方程式(10)~式(12)实现:
2.3 固壁边界处理方法
采用一种对Monaghan边界斥力法[20]进行改进的固壁边界处理方法叫几何平面边界。这种边界是一种虚拟的基于几何计算而确定的边界,采用这种边界处理方法能大量节省用来存储固壁边界粒子的内存,同时也提升了计算速度。这个几何平面由方程式(13)确定,
式(13)中,(a,b,c)是该平面的法向量,d来确定平面到坐标原点及坐标轴的距离。任意一个坐标为(x1,y1,z1)的水粒子到固壁平面的距离可由式(14)计算,计算出水粒子到固壁边界距离再根据Monaghan的斥力计算公式来计算固壁对水粒子的作用力。这种几何平面边界一般用来处理与坐标轴重合和平行的固壁边界。
3 模拟结果与分析
3.1 地下空间楼梯洪流基本情况
地下空间楼梯洪流数值计算模型是基于Ishigaki等[6]实验模型的物理尺寸进行建模,Ishigaki实验模型包含高水槽、阶梯以及从高水槽延伸到阶梯的地面平坦部分,模型尺度为13 m×2.5 m×6 m。台阶长度为0.3 m,高度为0.15 m,宽度为1 m,共20个台阶,高水槽中引入水泵往高水槽中持续加水。
为了突出GPU并行计算的计算效率,GPU并行计算的初始化条件与SPH串行模型[5]初始化条件稍有不同,具体数值建模及初始化如下:①采用基于虚拟几何平面的固壁处理计算区域外围四墙壁,布置固壁边界粒子来处理水槽的底部和楼梯台阶及高水槽到台阶的地面平坦部分,固壁边界粒子和水粒子的粒子间距均为0.037 5 m;②高水槽内初始水柱为4.5 m×2.5 m×2.0 m,水槽内水柱水平面比地面平坦部分高1.0 m,布置水粒子和固壁边界粒子共为496 563个,如图1中所示,蓝色区域为水粒子,绿色平面为固壁粒子;③初始时间步长为0.00003 s。GPU并行计算选用的SPH方法的数值处理技术与文献[18]中CPU串行模型的相同。
3.2 地下空间楼梯洪流模拟结果分析
利用改进的GPUSPH模型对地下空间楼梯流的计算结果如图2~图4中所示。高水槽内水柱发生坍塌后,沿地面平坦部分流向楼梯,从图3和图4中可以看出,采用与CPU串行计算相同的阶梯尺寸进行模拟计算时,也呈现了比较好的滑行流流态,这与Ishigaki实验结果和Gotoh的MPS方法的数值结果定性上比较一致。
同时,本文也对GPU并行计算的一些计算细节进行了统计,主要目的是与Gotoh的基于MPS粒子法的阶梯流模型和基于CPU串行计算的SPH模型进行对比,详细计算设置及计算结果见表1。SPH方法作为一种无网格粒子方法,虽然它比较适合处理具有大变形和移动交界面的问题,包括剧烈的阶梯流问题,但是作为一种使用粒子来离散计算区域和控制方程的无网格拉格朗日方法,它的计算量比较大,因为每个时间步都要计算每个粒子与其相邻近粒子的相互作用。因此,采用串行计算(CPU:Intel Core i7)时SPH的计算效率较低,然而SPH每个粒子的计算流程一致,采用GPU(NVIDIA Tesla C2075)加速后计算效率会显著提升。如表1中所示,基于GPU并行加速的SPH阶梯流模型计算效率非常高,对近50万粒子计算物理时间长度10 s只需要45 min,而同等计算尺度的CPU串行SPH模型却需要115 h,基于直接求解泊松方程的MPS模型计算更小的尺度,计算却更加耗时,需要296 h,计算效率最低,这也表明标准的SPH方法引入可压的状态方程来求解压力的重要性。
5 结论
本文采用SPH方法并借助GPU并行计算对地下空间阶梯流进行了模拟计算,并与基于MPS粒子法的串行计算结果和基于SPH方法的串行计算结果进行了对比分析,发现:
(1)通过引入可压状态方程求解压力场的串行SPH方法比串行MPS粒子法通过直接求解泊松方程来求解压力场的计算效率要高一倍多;
(2)经过并行化加速的SPH阶梯流模型计算效率非常高,比串行的SPH阶梯流模型和MPS阶梯流模型均要高出一个量级,相差不大的计算尺度情况下,串行的MPS阶梯流模型需要296 h,串行的SPH阶梯流模型需要115 h,而并行化的GPUSPH模型只需45 min。
(3)基于GPU并行化的SPH地下空间阶梯流模型不仅计算效率较高,而且能呈现较好的地下空间阶梯洪流流态,GPU并行化的SPH模型应用于较大尺度的城市街区洪水计算具有广阔的前景和可行性。
摘要：为了更好地模拟分析地下空间阶梯流动力学规律和提升数值计算效率,利用光滑粒子流体动力学方法(SPH),并借助GPU并行计算技术对地下空间阶梯流进行数值建模和计算。基于GPU并行化加速的SPH地下空间阶梯流模型的数值计算结果表明:该模型不仅计算获得较好的阶梯流流态,而且计算效率大幅度提高,比串行SPH模型和移动粒子半隐方法(MPS)模型快数百倍,为研究大尺度城市街区洪水提供了一个比较有前景的研究途径。
计算流体动力学仿真 篇7
为研究锂离子电池的热失控过程, FAA对可充电锂离子芯易燃性进行分析, 结果表明, 可充电的锂离子电池用很小的酒精火就可点燃[2]。王青松等讨论了锂离子电池热模型建立的两种途径, 通过热模型的建立, 指导锂离子电池的安全设计和管理[3]。吴凯等讨论了影响锂电池安全测试项目的因素, 并总结了改善锂离子电池安全性能的方法[4]。陈玉红等对锂离子电池爆炸机理进行分析[5]。然而这些研究主要集中于锂离子电池安全性能及对其进行安全改进, 对包装内锂离子电池发生热失控后的多米诺效应并无涉及。基于此, 为分析锂离子电池热失控的传播规律, 并研究其热失控传递的控制方法, 本文采用基于计算流体动力学 (CFD) 的数值模拟方法, 结合多米诺效应模型, 通过定量数据描述包装内锂离子电池热失控的多米诺效应, 分析锂离子电池间的热失控传播控制方法, 研究缓解锂离子电池热失控多米诺效应的控制方法。
1 锂离子电池热失控多米诺效应理论模型提出
1.1 锂离子电池热失控多米诺效应概述
多米诺效应常被用来分析事件之间一系列的连锁反应, 尤其是一个低概率事件可能引起的重大后果。多米诺效应被用于事故发生发展的理论研究也有相当长的时间。最广为人知的就是海因里希的多米诺模型, 也称为海因里希事故因果连锁理论, 是最早的顺序事故模型, 将事故解释为了一系列的离散事件按照一定顺序发生的结果[6]。但是海因里希的多米诺模型只是一个事故发生的概念模型, 并没有解释事故发生和危险出现的原因及其概率。
针对事故发生的多米诺效应的研究通常需要结合其他定量研究方法。例如, 在石油化工领域, 多米诺效应常与贝叶斯网络理论结合, 被应用于危险化学品储罐的事故发生可能性及其事故扩大后果的研究[7]。但是这种结合数学概率的方法要求已知研究对象的事故发生概率, 通常由事故数据经统计方法得到的经验模型而来。因此统计样本量及统计过程中样本之间的客观差异等其他因素都会对统计数据的结果可信性产生影响。采用仿真分析的方法进行事故发生多米诺效应的研究, 可以避免使用统计学方法时样本个体特性造成的偏差, 更接近事故发生的真实场景。
建立的锂离子电池热失控的多米诺效应理论模型见图1。将现实中多米诺骨牌倒下的过程抽象为电池热失控传递过程。将每个锂离子电池视为多米诺效应模型中的骨牌, 骨牌倒下代表相应电池发生热失控, 以锂离子电池热失控释放的能量为介质, 骨牌碰撞导致下一级骨牌代表锂离子电池热失控能量的传递。骨牌的连续倾倒的过程就是锂离子电池热失控的传递过程。
1.2 锂离子电池热失控多米诺效应理论模型建立
模型中假设初始热失控已经发生, 即第一级骨牌已经倒下。锂离子电池热失控的传递风险由风险集R确定。
式 (1) 中Ri代表i级骨牌倾倒导致i+1级骨牌发生热失控的风险。相邻两级骨牌发生多米诺效应的风险主要取决于两个因素, 骨牌的热失控释放能量E和骨牌间能量传递效率η。
式中, Ei代表第i级每个骨牌热失控释放能量的集合, ηi代表第i级骨牌热失控能量传递给i+1级每个骨牌的效率的集合, m为i级锂离子电池单元数量, n为i+1级锂离子电池单元数量。相应的, Eij代表第i级第j个骨牌热失控释放的能量。ηij代表第i级骨牌热失控能量传递给i+1级第j个骨牌的效率。
公式 (2) 中, Ri与Ei、ηi都成正比。因此控制锂离子电池热失控能量在某两级间传递的风险, 可以从降低热失控释放能量Ei和骨牌间能量传递效率ηi两个方面入手。根据锂离子电池的特点, 降低Ei的途径有:
(1) 降低单个骨牌热失控释放的能量。虽然理论上可以改进电极和电解液材料从本质上预防锂离子电池热失控的发生, 但是本质安全是不易实现的。减少单个电池的容量来减少Ei虽然可行, 却并不符合现实生产需要。
(2) 减少Ei集里元素的个数。从实际意义上讲就是减少参与放热的锂离子电池数量。因此在热失控传递过程中, 减少热失控热量影响范围内可能发生热失控的锂离子电池, 是降低多米诺效应的一种方法。是否可以通过减少某些位置极少量的锂离子电池有效的推迟热失控传递时间, 是一个值得深入研究的问题。
骨牌间能量传递效率ηi j可以用下面的方程表示:
式 (5) 中, dij表示i级骨牌与i+1级第j个骨牌的间距, hij表示i+1级第j个骨牌受撞击位置。d0代表热失控传递的临界距离, 超过此极限骨牌由于没有发生接触, 不发生传递现象。参数h0代表撞击有效范围, 对于锂离子电池就是它发生热失控的临界温度, 当到达下级锂离子电池的热量不足以使电池升至h0以上时传递效率为零。
结合公式 (4) 、式 (5) 和锂离子电池热失控发生的条件, 控制ηij有以下两种方式:
(1) 增加骨牌间距。第一, 增强电池之间的散热, 可以通过合理设计热管理系统加强散热, 或者在热失控初期被发现时采取有效的灭火散热手段实现。第二, 在骨牌之间设置能量屏障。降低热量传递的速度, 等效于增大了电池间距。增加隔热层, 使热失控释放的热量无法到达该级锂离子电池。是值得考虑的控制方案, 例如在电池组内增加电池间的隔热阻燃隔断层, 或者在包装内电池之间使用阻燃隔热的隔板。波音-787飞机APU使用的锂离子电池组就在电池之间增加了陶瓷的隔板以降低电池组整体火灾风险。
(2) 提高电池发生热失控的临界温度。同降低单个骨牌热失控能量一样, 这是不容易实现的。
对于上述控制方式可以通过实验研究或者数值仿真的方法来验证其有效性。对于锂离子电池热失控传递实验, 还存在电池平均温度测试的困难。而使用数值仿真可以方便的测量模型中任意位置的温度变化, 并且相比实验研究具有经济性好, 时间周期短, 控制单因素变量简便的优点。下面以空运包装内的18650型锂离子电池热失控多米诺效应控制仿真分析来验证理论模型的可信性。
2 锂离子电池热失控多米诺效应控制实例分析
2.1 空运锂离子电池包装介绍
由于锂离子电池的特殊火灾危险性[8,9], 航空运输锂离子电池有严格的包装要求, 除了瓦时额定值不能超过允许量以外, 还必须严格的包装在符合空运包装性能标准的内包装中, 才能装入外包装进行运输[10]。锂离子电池必须包装在完全封闭的内包装中, 并加以保护防止同一包装件内导电材料接触发生短路。
18650型锂离子电池是最常见的锂离子电池之一, 长度65 mm, 半径9 mm。其空运内包装情况如图2, 100个电池单元被装入一个专用纸盒中, 电池中间放置隔板以防止发电池倾倒和短路, 几何模型见图3。
分析锂离子电池热失控的多米诺效应理论模型知道, 对于运输中非工作状态的锂离子电池, 减少热失控传递途径上的少量电池和增强电池间能量屏障是降低热失控传递风险的有效方式。本研究以18650圆柱型锂离子电池在包装盒内的热失控传递为例, 结合仿真分析方法, 验证锂离子电池热失控多米诺效应模型得出的控制方式的可行性。
2.2 锂离子电池热失控多米诺效应的有限元模型
2.2.1 模型简化及网格划分
数值仿真采用的是商用CFD软件, 通过计算流体力学方法, 分析热失控在锂离子电池间的多米诺效应。由于锂离子电池包装中, 角落位置电池受到机械撞击和由于机械原因造成的内外部短路的可能性较大, 且研究热失控最初扩散更具有研究价值, 同时为了减少计算量, 选取包装内靠近角落的9个锂离子电池作为模拟分析的对象, 划分网格, 建立有限元模型 (见图4) 。为方便仿真分析结果的描述对电池编号 (见图5) , 其中编号为A的电池发生初始热失控, 中心位置电池编号为C, 距离电池A位置最远的电池编号为F。电池A、电池C和电池F位于该组电池的对角线位置, 关于对角线对称的其余电池分别编号为B1和B2, D1和D2, 以及E1和E2。
2.2.2 模拟中使用的材料参数及基本假设
模拟中考虑热传导和热辐射, 热对流在该模型中的贡献并不显著。使用的电池材料来源于文献数据及经验, 并不代表某一厂家或具体型号的电池 (见表1) 。锂离子电池热失控的放热曲线采用文献[11]中的锂离子电池随时间的典型放热曲线 (图6) 。模拟过程中假设:外部包装材料可以在计算时间内有效控制电池热失控放出的热量;电池在温度上升至180℃左右时达到热失控临界状态[12], 发生热失控;不考虑热失控放出可燃气体在包装内的扩散燃烧。
2.3 包装中锂离子电池热失控多米诺效应仿真结果及分析
2.3.1 减少中间骨牌对热失控传播的影响
监测电池B1~E2的平均温度, 达到180℃时按照热释放速率曲线加入热源, 观察电池到达热失控温度的时间。设定电池间使用具有阻燃能力的隔板, 隔板厚度为1 mm, 从电池A发生热失控放热开始, 262 s时电池F达到热失控温度, 此时的温度分布云图如图7所示。并且电池热失控多米诺效应骨牌的触发编号顺序是A>B1=B2>C>D1=D2>E1=E2>F。
依据锂离子电池热失控传递多米诺模型的分析, 减少某级骨牌数量可以降低锂离子电池热失控传播的风险。模拟中间位置电池C被阻燃材料替代的情况, 得到电池温度D1、E1、F的温升曲线, 与存在C电池的电池温度变化情况见图7。
从图8中可以看出, 使用阻燃材料替代电池C后, 电池D1热失控的间隔时间仅推迟了2 s, 但是电池E1和F的热失控开始时间有了较大幅度的延迟, 分别推迟了50 s、51 s。相对于电池C存在的情况, 电池F的热失控推迟了近20%, 可观的增加了电池热失控后的处置时间, 降低了热失控传递的风险。
电池D1热失控几乎没有被推迟的原因是热源B1、B2与电池D1相距较近, 即使消除C电池放热风险, B1、B2电池所在级别骨牌与D1所在级别骨牌的间距d依然在热失控传递的有效范围内, 因此热失控释放的热量依然可以直接对D1产生较大影响。
2.3.2 增加锂离子电池之间屏障的热阻力对热失控传播的影响
增加隔板厚度可以推迟下级电池的热失控时间, 可以推测, 当隔板增至一定程度, 将完全阻止热失控传递。考虑到实际包装条件, 在不改变包装内锂离子电池总数量的情况下, 允许的最大隔板厚度为1 mm, 因此模拟分析没有阻热隔板、0.5 mm隔板、1 mm隔板三种情况下C电池达到热失控临界温度的时间。仿真结果见图9。
由图9可知, 在没有隔板情况下C电池的达到热失控温度仅需78 s, 增加一个0.5 mm厚度的阻热隔板可以有效的推迟C电池发生热失控的时间多达98 s。继续增大隔板厚度至1 mm, 并没有大幅度的推迟热失控发生时间。因此增加阻热隔板对于控制锂离子电池的热失控传递有效, 但是限于当前包装布局和材料热阻的情况下, 增加的程度有限。
3 结论
本文建立了基于CFD仿真分析的锂离子电池热失控多米诺效应模型, 首先从模型入手定性剖析事故发生的可能性与控制方法, 进而通过仿真分析验证定性分析的准确性与热失控发展程度。分析结果显示两者结论具有一致性。结合空运锂离子电池包装的实际情况, 热失控传递控制的可行方法有:
(1) 减少运输包装中的极少量锂离子电池并使用阻燃材料代替其在包装中的位置, 就可以有效减缓热失控传递速度。
(2) 增加电池间的阻燃隔板, 可以明显的推迟热失控传递的时间。但是限于包装尺寸, 增大隔板厚度对推迟热失控传递的效果十分有限。
而推迟热失控传递时间对于及时扑救的有益效果以准确监测热失控发生为前提, 这有待今后进一步研究。
参考文献
[1] Federal Aviation Administration.Aviation cargo and passenger baggage incidents involving smoke, fire, extreme heat or explosion.http://www.faa.gov.2014-05-19
[2] Federal Aviation Administration.Flammability assessment of bulkpacked, rechargeable lithium-ion cells in transport category aircraft.DOT/FAA/AR-06/38, 2006
[3] Wang Qingsong, Ping Ping, Zhao Xuejuan, et al.Thermal runaway caused fire and explosion of lithium-ion battery.Journal of Power Sources, 2012;208:210—224
[4] 吴凯, 张耀, 曾毓群, 等.锂离子电池安全性能研究.化学进展, 2011;23 (2) :401—409Wu Kai, Zhang Yao, Zeng Yuqun, et al.Safety performance of lithium-ion battery.Progress in Chemistry, 2011;23 (2) :401—409
[5] 陈玉红, 唐致远, 卢星河, 等.锂离子电池爆炸机理研究.化学进展, 2006;18 (6) :823—831Chen Yuhong, Tang Zhiyuan, Lu Xinghe, et al.Research of explosion mechanism of lithium-ion battery.Progress in Chemistry, 2006;18 (6) :823—831
[6] Rausand M.Risk assessment thirty, methods, and applications.刘一骝译.北京:清华大学出版社, 2013:112—113Rausand M.Risk assessment thirty, methods, and applications.Liu Yiliu.ed.Beijing:Tsinghua University Press, 2013:112—113
[7] 张青松, 杨一楠, 郑薇.基于贝叶斯概率的油库池火灾多米诺效应分析.安全与环境学报, 2013; (03) :239—242Zhang Qingsong, Yang Yinan, Zheng Wei.Analysis on pool fire Domino effect in fuel tank district based on Bayes probability.Journal of Safety and Environment, 2013; (03) :239—242
[8] Federal Aviation Administration.Flammability assessment of lithiumion and lithium-ion polymer battery cells designed for aircraft power usage.DOT/FAA/AR-06/38, 2006
[9] Federal Aviation Administration.Risk in transporting lithium batteries in cargo by aircraft.SAFO 10017, 2010
[10] International Air Transport Association.Lithium cattery guidance document.http://www.faa.gov.2014
[11] Blake D.Battery heat release data.http://www.fire.tc.faa.gov.2013-05-13
计算流体力学在化学工程中的应用 篇8
一、计算流体力学在化学工程中的基本原理
计算流体力学是通过数值计算方法来求解化工中几何形状空间内的动量、热量、质量方程等流动主控方程, 从而发现化工领域中各种流体的流动现象和规律, 其主要以化学方程式中的动量守恒定律、能量守恒定律及质量守恒方程为基础。一般情况下, 计算流体力学的数值计算方法主要包括数值差分法、数值有限元法及数值有限体积法, 其也是一门多门学科交叉的科目, 计算流体力学不仅要掌握流体力学的知识, 也要掌握计算几何学和数值分析等学科知识, 其涉及面广。针对计算流体力学的真实模拟, 其主要目的是对流体流动进行预测, 以获得流体流动的信息, 从而有效控制化工领域中的流体流动。随着信息技术的发展, 市场上也出现了计算流体力学软件, 其具有对流场进行分析、计算、预测的功能, 计算流体力学软件操作简单, 界面直观形象, 有利于化学工程师对流体进行准确的计算。
二、计算流体力学CFD在化学工程领域中的应用研究
1. CFD在搅拌中的应用分析
由于搅拌槽内的流场的流动具有复杂性, 目前对搅拌槽等混合设备的设计和经验成分也采用理论计算的方式, 在化工领域中, 化工工业规模的反应器存在不均匀性等特点, 不均匀性随规模扩大而加重, 因此, 对搅拌槽内部流场进行研究是非常有必要的。虽然许多化学家对化工领域中的搅拌机槽内的流场进行了分析研究, 如Harvey等人采用二维模拟计算搅拌槽内流场的流体, 但随着技术的不断改革与发展, 计算流体力学的引进, 改变了以二维模拟的计算方式, 计算流体力学的方法不仅可以节约化工研究成本, 采用实验手段不能获得的数据, 计算流体力学方法也可以获得。Sun等人利用计算流体力学中的湍流模型计算了搅拌槽内的气液两相流动, 并且对其进行了三维模拟, 通过实验研究表明, 计算流体力学的数值模拟能有效的计算搅拌器上部的气体部分, 但是, CFD数值模拟也存在一定的缺陷, 不能有效模拟搅拌器底部区域。计算流体力学CFD与多普勒激光测速仪LDV有效结合, 可以对搅拌装置能更深入的研究, 其主要原因是多普勒激光测速仪测量的数据可以准确验证计算流体力学CFD计算的结果, 同时多普勒激光测速仪测定特定点的速度也可以作为计算流体力学计算的参考条件。
2. CFD在化学工程换热器中的应用分析
换热器是化学工程中使用最多的设备, 通过计算流体力学的计算方式, 不仅可以精确、详细的测量换热设备内流场的流动, 也可以预测换热器的性能, 经济可靠的换热器对化工工业具有重要作用。对于化工中的管壳式换热设备, 其内部的几何形状设备结构复杂, 利用计算流体力学模拟管壳式换热设备的壳侧流场, 进而充分了解管壳式换热设备的壳侧在瞬间变化中的温度场、速度场, CFD的应用有利于分析研究换热器的基本原理和结构构造。
3. CFD在化学反应工程中的应用研究
针对高温的太阳能化学反应器, 温度、压力场和计算速度等数据指标等可以从CFD模拟优化中获取, 包括粒子的运动轨迹, 但当温度大于1500k和辐射强度超过3000KW/m2的状态时, 这些数据是无法获取的。学者采用有限元法的计算机流体力学软件的计算方式, 对氧化燃烧化学反应中的催化剂进行模拟, 最后, 分析实验结果。由于聚乙烯是化学反应中使用最多的聚合物, 其化学反应也需要在高温、高压的环境下通过自由基聚合的方式生成, 但是聚乙烯的生成需要利用绝热高压釜反应器, 学者为了更深入地研究绝热高压釜反应器的内部构造和化学反应原理, 对反应器中的旋转浆叶、湍流模型等采用CFD中的FLUENT软件进行三维模型, 经过计算, 可以得出不同环境下的反应器的流线、反应器内部的浓度梯度及温度梯度。通过CFD软件预测反应器的速度、温度及压力场, 可以更进一步理解化学反应工程中的聚合过程, 详细、准确的数据可以优化化学反应中的操作参数。
结束语:
计算流体力学在化学工程被得到广泛的应用, 大多数化学工程师及专家都纷纷采用计算流体力学CFD研究化学领域中内部流场的流动流体。因此, 加强计算流体力学软件的研究, 使计算流体力学应用在化学工程中的各个领域。
参考文献
[1]于艳, 樊耀波, 徐国良, 王琦, 杨文静.计算流体力学在膜技术及膜生物反应器研究中的应用[J].膜科学与技术, 2011, 01:105-112.
[2]刘国钰.计算流体力学在优化清水池水力中的应用[J].产业与科技论坛, 2012, 20:68-69.
[3]于艳, 樊耀波, 徐国良, 王琦, 杨文静.计算流体力学在膜技术及膜生物反应器研究中的应用[A].中国膜工业协会、同济大学环境科学与工程学院、中国城镇供水排水协会、上海市政工程设计研究总院、哈尔滨工业大学.2010年膜法市政水处理技术研讨会论文集[C].中国膜工业协会、同济大学环境科学与工程学院、中国城镇供水排水协会、上海市政工程设计研究总院、哈尔滨工业大学:, 2010:9.