星期三, 一月 03, 2007

推荐几本学ANSYS的好书

这里要推荐四本学ANSYS的好书。这当中前两本是我经常翻看的,另外两本是我今天在一个ANSYS用户那儿看到,觉得很好,对于一些ANSYS高级用户可能会有用。这些书都是英文的,我给出了Google book search和Amazon的链接。这些书的价格都是很贵的,如果要买的话,需要有老板赞助。

当然最好的ANSYS学习参考书是ANSYS的在线帮助。

1. Finite Element Analysis: Theory and Applications with ANSYS, Second Edition (有限元分析:理论及ANSYS应用,第二版)
作者:Saeed Moaveni

Google book search链接:无,没有收录
Amazon链接:http://www.amazon.com/Finite-Element-Analysis-Theory-Applications/dp/0131112023

简介:这本书从最基本的有限元分析方法讲起,涵盖了基本的有限元单元的推导(truss,beam,plane结构单元和1d,2的传热单元),并且还简单介绍了模态分析。这本书适合初学者(包括没有任何有限元分析基础的人)。书里面的例子简单易懂,推导循序渐进,并且理论推导结合简单的手工计算实例和ANSYS计算实例。手工计算实例和ANSYS计算实例的讲解非常详尽。是一本学习有限元基础知识和ANSYS基本操作的好书。


2. An Introduction to Computational Fluid Dynamics: The Finite Volume Method Approach (计算流体力学简介:有限体积法)
作者:H. Versteeg, W. Malalasekra

Google book search链接:无,没有收录
Amazon链接:http://www.amazon.com/Introduction-Computational-Fluid-Dynamics-Approach/dp/0582218845

简介:这本书同样是一本入门书籍,非常适合于初学者。和上面一本一样,这本书对基本概念的阐述非常详尽和清楚,公式推导过程清晰易懂。该书简介CFD的基本思想,和湍流模拟的基本概念。该书对有限体积法的离散方法的作了详细介绍。事实上,认真研读本书足以编出基于结构化网格的简单CFD程序。

题外话:这本书其实不是ANSYS参考书,但是由于上面的两本书分别是我学FEA和CFD的启蒙书,加上他们的编排比较适合初学者,所以一起放在一起。对于通常的FEA和CFD分析,掌握这些基本知识可能就足够了。另外,对于CFD,陶文铨教授的《数值传热学》我也认真看过。可惜一直没有能买一本。《数值传热学》是我看到的国内教授写的最好最有分量的一本书。听说陶老师刚刚被评为中科院院士,恭喜。另外,这部CFD书的作者H. Versteeg, W. Malalasekra中的一个是我一个同事的毕业论文导师。据他说,这两个人一个好,一个一般,哈哈。

3. The Finite Element Method and Applications in Engineering Using ANSYS (有限元方法及ANSYS的工程应用)

作者:Erdogan Madenci, Ibrahim Guven

Google book search链接:无,没有收录
Amazon链接:http://www.amazon.com/Finite-Element-Method-Applications-Engineering/dp/0387282890

简介:上面我说过学习ANSYS最好的参考书是ANSYS的在线帮助。而如果翻阅一下这本书的目录,你会发现这本书的结构和ANSYS的在线帮助有几分相似。通常的有限元书籍,可能会花很多篇幅讨论有限元方程的建立,但是对于工程中常用的非线性材料模型,接触分析等很少有分析。因此,这本书适合对ANSYS有一定了解,向知道一些进阶功能的朋友。


4. Fatigue Life Prediction of Solder Joints in Electronic Packages with ANSYS (ANSYS在电子封装焊点疲劳分析中的应用)
作者:Erdogan Madenci, Ibrahim Guven, Bahattin Kilic

Google book search链接:http://books.google.com/books?vid=ISBN1402073305&id=E7HOfOQsgpEC
Amazon链接:http://www.amazon.com/Prediction-Electronic-Packages-International-Engineering/dp/1402073305

这本书对于电子组件焊点在温度循环(thermal cycling)条件下的疲劳寿命的分析非常有用。焊点疲劳分析是一种常见的电子封装有限元模拟。该书不仅回顾了焊点疲劳分析的理论,而且提供了非常用户友好的宏。通常做温度循环分析,需要做submodeling。而且电子组件的结构非常复杂,需要很长的时间(几个星期)才能生成高质量的六面体网格。该书提供的宏和图形界面模块,很好的解决这样的问题。用户只需输入焊点以及其他部件的大小尺寸,系统会自动生成六面体网格,并且进行温度循环分析。并且有专门的单元类型切换,submodeling,材料定义,后处理模块,非常强大和人性化。

这本书的价格虽然很贵(186美元),但是该书所附的宏,绝对是物超所值。

星期一, 一月 01, 2007

ANSYS/Fluent年会论文下载

最近安世亚太飞昂软件都相继把他们的年会论文放在网上,供免费下载。

2006 年ANSYS用户大会暨公司十年大庆

官方网站:http://www.ansys.com.cn/10th/

论文下载:http://www.ansys.com.cn/10th/paper.html

2006中国FLUENT用户大会

官方网站:http://www.fluentchina.com/cfd_index.html

论文下载:http://www.fluentchina.com/news_news_20061220.html

另外一个消息是,ANSYS11.0到目前还没有发布,看来台湾的虎门科技是摆了乌龙了。最近几日,在打开的workbench时,左边的新闻栏会显示,ANSYS和IE7的兼容性问题,当中隐约提到ANSYS11.0会在二月发布。这也算是正常的,现在是圣诞假期,怎么会发布新版呢?

FDS、Smokeview、PyroSim入门教程

FDS是Fire Dynamics Simulator(火灾动力学模拟器)的缩写。官方网站是http://www.fire.nist.gov/fds。FDS是一个专注于火灾模拟的CFD求解器,而Smokeview则是它的后处理器。FDS本身是没有图形化前处理器的,也就是用户要先编辑命令定义要求解的问题。FDS是由美国NIST(National Institute of Stardard and Technology,美国国家标准与技术研究院)开发的免费(对的,免费,因此你可以下载安装)CFD软件。Pyrosim是我在网上找到的(Google里搜索关键字FDS fire,右边会显示赞助商链接)。它是第三方开发的图形化FDS前处理器。PyroSim不是免费软件,但是提供30天的试用。对于学习者来说,30天应该是足够了。

想在网上找一些关于FDS的中文资料,但是系统的介绍很少。这个博客(人间四月芳菲尽)对FDS有零星的介绍。可惜的是,有一段时间没有更新了。在sogou里面搜索Pyrosim,只有屈指可数的几个结果。看来这个软件的应用还不是很广泛。

FDS的下载地址是:http://www.fire.nist.gov/fds/download/fdsall_407a.exe。在官方网站上,还有其它相关信息,比较有用的是FDS/Smokeview example,当中有超过150个算例,可以让初学者了解一下软件的功能。其他的文档在官方网站或者安装目录下都可以找到。其中,User's gudie(用户手册)是介绍软件的基本使用,而Technical reference guide(技术参考手册)是介绍求解器的模型及背景理论,需要一定的CFD基础以及相当长的时间消化。

对于FDS的input file(输入文件),这里就不多做解释了,一是因为本人并不熟悉这些功能,而是因为PyroSim提供了图形化前处理功能,把用户从枯燥的命令行中解放了出来。Pyrosim的官方网站是http://www.pyrosim.com/。用户提交了基本注册信息后,会收到一封含有下载地址的email。安装文件大概有45M。实际上,该安装文件不光有Pyrosim主程序和帮助,还有FDS/smokeview的主程序及文档。也就是说,如果你决定使用PyroSim,是不用单独下载FDS的。

下面我们要来做两个模拟,一个是Pyrosim Example Guide里面的第一个例子(Chapter 1. Air Movement),然后修改这个例子,做一个隧道火灾的模拟。

1. 运行一个FDS例子

1.1 启动PyroSim(并打开Example Guide, Chapter 1)。具体就不多说了,开始>程序>…………

下面是PyroSim启动时的截图,上面有中文“火”的写法,但是看起来似乎是源自日文汉字(Kanji)的火。但是开发PyroSim的Thunderhead engineering是一家位于美国的公司。

1.2 建立网格

选择Model > Edit Grid ... 在跳出的面板上点New。把网格的边界设置为Min X = 0,Max X = 10,Min Y = 0,Max Y = 10,Min Z = 0,Max Z = 10。同时把网格数量设为X cells = 20,Y cells = 20,Z cells = 20。如下图所示:

在FDS里面,基本形状只能是长方体的(尽管通过长方体的组合,能模拟复杂的几何形状)。网格也只能是正交网格(点确定后可以看到)。这个例子是模拟风洞内的燃烧,风洞尺寸为10米*10米*10米。每边划分成20个网格,总共有8000个网格。在FDS中,由于求解器的特殊要求,每边划分的单元数通常要求是2,3或5的倍数。

点击OK后,网格就生成了,如下图所示:

1.3 定义粒子

这里的粒子并不是通常CFD软件里的拉格朗日粒子,只是为了后处理方便而定义的示踪粒子。选择Model > Edit Particles ...,然后点New。不要选Particle Have Mass,这样粒子没有质量,仅是示踪粒子。选中Color Particles During Animation,使用默认的颜色,红色。这样在结果中会有红色的示踪粒子。

1.4 建立面

实际上,在建立网格的时候,我们已经得到了计算模型。这里的面是用来定义边界条件。要注意的是,这里仅仅定义边界条件,而没有给模型中的面指定边界条件。也就是说只定义边界上的物理条件,但并没有和模型中的边界联系起来。

定义入口边界。Model > Edit Surface Properties ...,点New。Surface name是BLOW,使用INERT作为模板。选择surface type为 "Fan/Wind”。然后下面会有很多选项。系统定义的表面类型都有明确的物理意义,还是比较好懂的。Fan/Wind里可以定义入口空气的温度和速度。为了便于区分,把Color选成蓝色。Air Temperature = 20 C,Specify Normal Velocity = -1.0 m/s。-1.0 m/s代表气流速度为1m/s,方向为进入计算域(+1.0m/s为出)。第三个标签下,选中Emit Particles,Particle Type = PART。最后点OK。各步骤的截图如下:

在这里例子里面,只需要定义BLOW这一个面,其余的面可以用系统定义的来表示。

1.5 定义边界条件

这里是真正指定表面的边界条件。

定义入口。选择Model > New Vent ...,在Specification标签下,Description = Vent Blow,Type = BLOW,Lies in the plane X = 0.0,Min Y = 3,Max Y = 7,Min Z = 3,Max Z = 7。在X=0平面上,有4*4的范围为BLOW的边界条件。点OK,离开界面。具体操作和操作后的结果为:

定义出口。步骤和上面的基本一样,Model > New Vent ... > New。Description = Vent Open,Type = Open,Lies in Plane X = 10,Min Y = 3, Max Y = 7,Min Z =3,Min Z = 7。在和入口相对的面上,有一个4*4的通风口。

1.6 定义切面

FDS使用LES模型模拟湍流,会产生大量的瞬时数据。FDS里面需要定义切面,只有切面数据会保存。

Output > Slices ...。XYZ plane = Y,Plane value = 5,Gas Phase Quantity = Velocity,Use Vector = No。OK。在Y=5平面上,保存速度值,但是不保存矢量。

1.7 设定模拟参数

FDS > Simulation Parameters ...。Simulation Title = Wind Tunnel,Specify Duration = 60 s,Initial Time Step = 0.15 s。数值模拟总时长为60秒,初始时间步长为0.15秒(后面的时间步长有系统自动决定)。OK。

1.8 运行FDS

FDS > Run FDS ... 先要保存一个*.data文件(自行命名)。然后FDS开始求解,求解过程如下图所示:

1.9 检查结果

FDS算完60秒后,会自动弹出SmokeView窗口。 在上面点右键Load/Unload > Slice File > VELOCITY > *Y=5.0。会显示前面定义的切面上的速度云图。

在上面点右键Load/Unload > Particle file > *SMOKE/WATER。会显示粒子轨迹。

至此,我们已经完成了FDS里面的第一个练习。下面可以修改上面的例子,做一个隧道内火灾的模拟。

2. 隧道火灾模拟

隧道内的火灾模拟是一个很常见的FDS应用。如果在隧道内发生火灾,如果通风系统没有工作,燃烧产生的烟会向两边对称扩散,引起能见度的下降,和对隧道内人员的健康威胁。通常,需要在隧道入口装风扇,把烟吹向出口,这样烟雾就不会向入口扩散。烟雾向上游的扩散成为是逆流(back layer flow)。入口风速大,逆流长度就小,或者没有逆流;入口风速小,逆流长度就长。逆流长度刚好为零的入口风速成为是临界速度(critical velocity)。FDS广泛使用于逆流的研究。

这里要做的是一个60*16*6的隧道。入口速度是3m/s,入口温度是25C。在底部有一块8*8的燃烧区域,火灾大小为100MW。为了简化问题没有模拟过程,而是直接把100兆瓦的热量直接分布在64平方米的面积上,折合1562.5kw/m2。

2.1 继续前面的模型

下面要修改前面建立的模型。如果模型没有打开,可以打开保存的*.psm或者*.data文件。

2.2 建立网格

在左边的树形目录上,双击Grids > GRID。在弹出的面板上,修改Max X = 60,Max Y = 16,Max Z = 6,X cells = 120,Y cells = 32,Z cells = 12。OK。点OK后,可能模型不在屏幕中央,可以点击工具栏上右数倒数第二个图标(Reset View to All Visible Objects)重置显示。步骤和结果如下所示:

2.3 建立面

这里需要建立两个面。第一个面是入口,修改前面的BLOW条件;第二个面是FIRE,需要创建。
在左边双击Surfaces > BLOW。Properties > Air Temperature = 25 C。Air Flow > Specify Normal Velocity = -3 m/s。

Model > Edit Surface Properties ... > New...,Surface Name = FIRE,OK。Surface Type = Non-Flammable Solid。Boundary Conditions > Boundary Types = Fixed Heat Flux,Heat Flux = 1562.5,OK。

2.4 定义边界条件

按照上面的描述,有三个边界条件需要定义。

双击Model > Vent Blow。把整个X=0平面设成是入口。Min Y = 0,Max Y = 16,Min Z = 0,Max Z = 6。OK。

双击Model > Vent Open。把整个X=60平面设成是出口。Lies in the plane X = 60,Min Y = 0,Max Y = 16,Min Z = 0,Max Z = 6。OK。

新建一个Fire Region边界条件,在主菜单上选择:Model > New Vent。Description = Vent Fire,Type = FIRE,Lies in the plane Z = 0,Min X = 26,Max X = 34,Min Y = 4,Max Y = 12。OK。

过程及结果如下:

2.5 定义切面

在Y=8平面定义一个切面。Output > Slices ...。把Plane Value改成8。在Y=8平面定义一个新的切面,来显示温度。因为我们没有模拟燃烧,根据温度场可以大概知道烟雾的分布。XYZ Plane = Y,Plane Value = 8,Gas Phase Quantity = TEMPERATURE,Use Vector? = No。

2.6 设定模拟参数

模拟的总时长为60秒。FDS > Simulation Parameters ...。Simulation Title = Tunnel Fire,Specify Duration = 60,Initial Time Step = 0.1。OK。

2.7 运行FDS

FDS > Run FDS ...。先保存一个合适的*.data文件。然后求解器会启动。求解过程大约为15分钟(在我的电脑上)。求解结束后,会自动跳出SmokeView。

1.9 检查结果

在自动弹出SmokeView窗口上面点右键Load/Unload > Slice File > TEMPERATURE > *Y=8.0。会显示前面定义的切面上的温度云图。

1.10 分析

在上面的温度云图中,我们可以看到,在入口速度为3m/s的时候,隧道内有非常显著的回流。也就是说3m/s的入口速度,不足以阻止烟气向上游扩散。有兴趣的朋友可以增加入口速度,看看这个问题的临界速度是多少。我试出的临街速度是4.5m/s左右。

另外一个问题是这个结果准确吗?理论上来说LES和k-e相比,有很大的优越性,但是由于FDS使用是有限差分的方法,而且没有检查收敛性(explicit),通常情况下只能作为参考。比如这个例子,有兴趣的朋友,可以用Fluent/CFX建模算一下,Fluent/CFX给出的临界速度应该在2.7m/s左右。差别还是相当大的。另外根据我的测试(和一组条件不同的实验数据的比较),基于k-e的Fluent/CFX结果通常会稍稍低估临街速度。但是总体而言,还是Fluent/CFX的结果更加可信。值得提醒的是,用Fluent/CFX的时候,一定要记得加上浮力修正。

星期一, 十二月 25, 2006

什么是64位计算?

64位(64 bit)高性能计算是当前一个很热的话题,很多机构已经配置了、或将要升级至64位计算。但是实际上“64位”这个概念对我来说还是很模糊的。当然还有容易令人混淆的字(word),字节(byte),位(bit)的概念。

下面是我从互联网上找到的关于64位计算的一些资料:

1. Wikipedia: 64-bit

2. Microsoft: Overview of Windows XP Professional x64 Edition

3. IT168: 64位技术

对于64位技术,有各种不同的定义,Wikepedia的定义为:如果没有特别的描述,通常“64位” 计算机架构是指CPU整型寄存器有64比特的数据宽度。“64位”架构支持长度为64比特的整型数据。(Without further qualification, however, "64-bit" computer architecture generally has integer registers that are 64 bits wide, which allows it to support (both internally and externally) 64-bit "chunks" of integer data.)

64位计算主要优点是它能够支持更大的内存。由于内存地址是由整型数表示的。更大的整型数便意味着更大的内存支持。理论上,32位整型寄存器支持2^32个内存地址,即4GB,而64位整型寄存器支持2^64个内存地址,即16TB。实际上,由于其他的限制,Windows 32位操作系统对单个进程支持2GB内存,目前使用的Windows x64操作系统支持128GB内存。

我们并不能简单的认为64位的计算性能是32位的两倍。通常64位的性能会稍稍高于32位,但同时占用的内存也会稍稍多一些。笼统地说,64位计算极大地扩展了内存地址,但是和计算速度,内存使用(比如浮点数长度,单精度/双精度)并没有显著关系。

另外,要使用64位技术必要要有64位处理器、64位操作系统和64位软件。三者缺一不可。而现在64位应用软件相当匮乏,虽然最新版本的CAE软件中很多都实现了对64位计算的支持。

对64位操作系统而言,1 word = 64 bit, 1 byte = 8 bit。因此1 word = 8 byte。1GB内存可以换算成1024 * 128个字:1 Giga byte = 1024 Mega byte = 1024 * 1024 byte =1024 * 128 word。

对32位操作系统而言,1 word = 4 byte,1 GB = 1024 * 256 word。

圣诞快乐

圣诞快乐!并附上一个朋友用ANSYS Workbench做的图片(一个简单但有趣的应力分析)。

星期六, 十二月 23, 2006

CD-adapco发布STAR-CD 4.0.2

新闻链接

从新闻稿看,增加的新功能并不是太多,比如前处理功能的更新(尤其是多面体网格功能)。有趣的是,Fluent 6.3和CFX的多面体网格都是基于常用的四面体或六面体网格,在求解器内部自动重新组合。也就是说不需要专门的网格划分软件。但是CD-adapco的多面体网格解决方案似乎要“笨”一些,需要专门的多面体网格划分软件。这样做的优点是什么呢?
当中重要(或者是最主要)的新功能应该是多物理场(mutli-physics)功能。虽然多物理场已经不是一个新概念(ANSYS Multiphysics),但是这次好像是第一次由CFD开发商提出来。 STAR-CD 4.0.2将可以进行固体应力分析和固体融化、固化分析。不知道STAR-CD的应力分析模块现在有哪些功能,是基于有限体积法还是基于有限元方法。

有限元法、有限差分法和有限体积法的区别

昨天在网上看到一篇关于有限差分,有限元和有限体积法区别的文章。现在贴在下面,也许对有兴趣的朋友有帮助(我对文章作了一些编辑)。

1. 有限差分法(Finite difference method)

有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域。有限差分法以泰勒(Taylor)级数展开等方法,把控制方程中的导数用网格节点上的函数值的差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。对于有限差分格式,从格式的精度来划分,有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑,可分为中心格式和逆风格式。考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式,主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分方法主要适用于结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。

构造差分的方法有多种形式,目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等, 其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合,可以组合成不同的差分计算格式。

2. 有限元法(Finite element method)

有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内,选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式 ,借助于变分原理或加权余量法,将微分方程离散求解。采用不同的权函数和插值函数形式,便构成不同的有限元方法。有限元方法最早应用于结构力学,后来随着计算机的发展慢慢用于流体力学的数值模拟。

在有限元方法中,把计算域离散剖分为有限个互不重叠且相互连接的单元,在每个单元内选择基函数,用单元基函数的线形组合来逼近单元中的真解,整个计算域上总体的基函数可以看为由每个单元基函数组成的,则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中,常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同,有限元方法也分为多种计算格式。从权函数的选择来说,有配置法、矩量法、最小二乘法和伽辽金法,从计算单元网格的形状来划分,有三角形网格、四边形网格和多边形网格,从插值函数的精度来划分,又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。对于权函数,伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ;最小二乘法是令权函数等于余量本身,而内积的极小值则为对代求系数的平方误差最小;在配置法中,先在计算域 内选取N个配置点 。令近似解在选定的N个配置点上严格满足微分方程,即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成,但也有采用三角函数或指数函数组成的乘积表示,但最常用的多项式插值函数。有限元插值函数分为两大类,一类只要求插值多项式本身在插值点取已知值,称为拉格朗日(Lagrange)多项式插值;另一种不仅要求插值多项式本身,还要求它的导数值在插值点取已知值,称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标,有对称和不对称等。常采用的无因次坐标是一种局部坐标系,它的定义取决于单元的几何形状,一维看作长度比,二维看作面积比,三维看作体积比。在二维有限元中,三角形单元应用的最早,近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元,常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。

对于有限元方法,其基本思路和解题步骤可归纳为

(1)建立积分方程,根据变分原理或方程余量与权函数正交化原理,建立与微分方程初边值问题等价的积分表达式,这是有限元法的出发点。
(2)区域单元剖分,根据求解区域的形状及实际问题的物理特点,将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作,这部分工作量比较大,除了给计算单元和节点进行编号和确定相互之间的关系之外,还要表示节点的位置坐标,同时还需要列出自然边界和本质边界的节点序号和相应的边界值。
(3)确定单元基函数,根据单元中节点数目及对近似解精度的要求,选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的,由于各单元 具有规则的几何形状,在选取基函数时可遵循一定的法则。
(4)单元分析:将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近;再将近似函数代入积分方程,并对单元区域进行积分,可获得含有待定系数(即单元中各节点 的参数值)的代数方程组,称为单元有限元方程。
(5)总体合成:在得出单元有限元方程之后,将区域中所有单元有限元方程按一定法则进 行累加,形成总体有限元方程。
(6)边界条件的处理:一般边界条件有三种形式,分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件,一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件,需按一定法 则对总体有限元方程进行修正满足。
(7)解有限元方程:根据边界条件修正的总体有限元方程组,是含所有待定未知量的封闭方程组,采用适当的数值计算方法求解,可求得各节点的函数值。

3. 有限体积法(Finite volume method)

有限体积法(Finite Volume Method)又称为控制体积法。其基本思路是:将计算区域划分为一系列不重复的控制体积,并使每个网格点周围有一个控制体积;将待解的微分方程对每一个控制体积积分,便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分,必须假定值在网格点之间的变化规律,即假设值的分段的分布的分布剖面。从积分区域的选取方法看来,有限体积法属于加权剩余法中的子区域法;从未知解的近似方法看来,有限体积法属于采用局部近似的离散方法。简言之,子区域法属于有限体积发的基本方法。

有限体积法的基本思路易于理解,并能得出直接的物理解释。离散方程的物理意义,就是因变量在有限大小的控制体积中的守恒原理,如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。 限体积法得出的离散方程,要求因变量的积分守恒对任意一组控制体积都得到满足,对整个计算区域,自然也得到满足。这是有限体积法吸引人的优点。有一些离散方法,例如有限差分法,仅当网格极其细密时,离散方程才满足积分守恒;而有限体积法即使在粗网格情况下,也显示出准确的积分守恒。就离散方法而言,有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律(既插值函数),并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值,这与有限差分法相类似;但有限体积法在寻求控制体积的积分时,必须假定值在网格点之间的分布,这又与有限单元法相类似。在有限体积法中,插值函数只用于计算控制体积的积分,得出离散方程之后,便可忘掉插值函数;如果需要的话,可以对微分方程中不同的项采取不同的插值函数。

星期四, 十二月 07, 2006

ANSYS 11.0终于即将发布

最近一两个月实在是太忙了,发文章很不正常,真是太内疚了。:(

这几天正在参加一个conference,晚上回来查看网页,无意中看到台湾虎门科技的网页(虎门科技是ANSYS在台湾的代理商)。居然得到了ANSYS即将发布的消息。台湾的发布会会在12月27号进行,也就是说ANSYS 11.0会在2-3周内发布。

星期一, 十二月 04, 2006

Fluent 6.3已经发布

昨天早上打开ANSYS的主页,发现上面的新闻稿说Fluent 6.3已经发布了。下面一个问题是ANSYS 11.0什么时候发布呢?现在ANSYS 11.0和10.0之间至少已经间隔16个月了。

Fluent6.3发布的相关新闻链接在这里这里是关于Fluent 6.3新的特性的概述。

1. 网格、数值及并行计算:新的基于压力的耦合求解器,新的建模检查工具,新增加到多面体网格的支持,增加对64位windows的支持以及其它并行计算能力的提高。

F1赛车多面体网格

2. 传热、相变和辐射:面对面辐射模型应用范围和效率的提高,用户可以耦合离散坐标辐射模型和能量方程,支持多组分实际气体。

玻璃熔炉模拟

3.动网格:动网格可以应该在稳态问题中,动网格的定义更加直接和方便,滑移网格支持一边多个网格。

4.反应流:新的低速化学反应和微混合模型,更多的化学物质和反应可以使用无预混合或部分预混合模型模拟,更多模拟SOx、NOx生成的模型。

5.湍流和声学:更灵活的DES模型,LES模拟更多的时间统计结果,提高了RSM模型,用户可以自定义壁面函数。

6.多相流:瞬态多相流模拟精度得到了提高,欧拉、混合和VOF模型和壳传热模型兼容,扩展了欧拉模型的应用范围,在使用VOF模型时,可以使用新的界面跟踪模型,可用用UDF设定自由液面接触角,对离散相模型可以模拟多组分例子的气化,在使用动网格时,粒子可以从动网格表面射入。

7.后处理:一些新的报告功能,增强了和第三方后处理软件的协作,Fluent并行计算结果可以输出到tecplot。

8.自定义工具:扩展了用户定义标量的应用范围,更多的UDF链接功能,对扩展模块(add-on modules)的改进。

星期日, 十一月 12, 2006

ANSYS计算节点内力的宏

今天根据ansys.net上的一个宏稍微修改了一下,做了一个计算节点内力的宏,贴在下面。这个宏的原始作者是Shen-Yeh Chen。这里是原始描述原始宏

这个宏可以计算指定节点上的内力。通常情况下,ANSYS只计算被限定自由度的节点的反应力。如果模型里有两个材料,并且这两个材料用连续网格模拟(而不是接触),有些时候我们想知道这两个材料间的相互作用力,而在ANSYS里面并没有直接的命令。使用这个宏可以做到这一点。

这个宏的使用法方法是:0.当然是把这个宏命名并放在相应的目录下;1.选择要计算内力的节点;2.选择这些节点一边的单元,因为这些节点代表一个截面,而这个截面会有两边有单元。选择不同的边,内力的方向是不同,如果两边都选的话,内力就会平衡为零了。因此必须要先选择一边,但是不必要只选和节点相邻的一层,因为宏里面已经有命令去选择和节点相连的单元;3.运行宏;4.在工作目录下,会有一个x_nfor2f.ans文件,这个文件的最后一行,就是在x,y,z三个方向的内力分量了。

这个宏必须要在通用后处理器中运行。

! ========================================

! Shen-Yeh Chen

! 05/01/1998

! 12/10/2001

! =========================================

!

! This macro will transfer the element force on the selected

! node into cross-section force on each node, and write

! all the force on each node into the "F" command. Note that

! the macro go through each SELECTED node, sum up all SELECTED

! elements that are attached to the SELECTED node. Then the

! force on each nodes is added a negative sign and printed

! into the 'F' command.

! ---------------------------------------------------------

! Syntax :

! x_nfor2f

!

! Input :

! None

!

! Output :

! file : x_nfor2f.ans : nodal force for the section load

! x_nfor2f.cdb : nodal force & nodal coordinates/ID in

! CDB format

! nodal components : x_nfor2f : defined in x_nfor2f.cdb

! X_NFOR2F_FX : summation of FX

! X_NFOR2F_FY : summation of FY

! X_NFOR2F_FZ : summation of FZ

! Temporary components

! X_NFOR2F

! X_NFOR2F_E

!

! Temporary components

! x_nfor2f.db

!

! Temporary variables

! AR30,AR31,AR32,AR33

!

! External Programs called

! None

!

!

! ========================================

/POST1

*GET,ZNNODE,NODE,,NUM,MAX

CM,X_NFOR2F,NODE

CM,X_NFOR2F_E,ELEM

ESLN,R

!

AR30=0

*CFOPEN,x_nfor2f,ans

*VWRITE

('/PREP7')

X_FOR2F_FX=0.0

X_FOR2F_FY=0.0

X_FOR2F_FZ=0.0

!

! LOOP THROUGH ALL SELECTED NODES

:LINE01

AR30=NDNEXT(AR30)

*IF,AR30,EQ,0,THEN

*GO,:LINE02

*ENDIF

NSEL,R,NODE,,AR30

FSUM

!

!

! GET THE FORCE

*GET,AR31,FSUM,0,ITEM,FX

*GET,AR32,FSUM,0,ITEM,FY

*GET,AR33,FSUM,0,ITEM,FZ

!

AR31=-AR31

AR32=-AR32

AR33=-AR33

X_FOR2F_FX=X_FOR2F_FX+AR31

X_FOR2F_FY=X_FOR2F_FY+AR32

X_FOR2F_FZ=X_FOR2F_FZ+AR33

!

! WRITE TO FILE

*VWRITE,AR30,AR31

('F,',F6.0,',FX,',E12.5)

*VWRITE,AR30,AR32

('F,',F6.0,',FY,',E12.5)

*VWRITE,AR30,AR33

('F,',F6.0,',FZ,',E12.5)

!

! SELECT THE NODES

CMSEL,,X_NFOR2F

*GO,:LINE01

!

!

! SAVE THE STATE TO x_nfor2f.db

:LINE02

*VWRITE,X_FOR2F_FX,X_FOR2F_FY,X_FOR2F_FZ

('FX = ',E12.5,', FY = ',E12.5,', FZ = ',E12.5)

*CFCLOS


星期一, 十一月 06, 2006

STAR-CD和STAR-CCM+ :冠军的选择

最近CD-Adapco网站上登出了一则新闻。新闻的题目是:STAR-CD and STAR-CCM+: Still the choice of champions (STAR-CD和STAR-CCM+:冠军的选择)。这则新闻被很多网站转载。新闻说的是2006赛季F1冠军车队雷诺车队使用的是CD-Adapco提供的CFD软件。我每天上班的路上看的一份小报几乎每天都要介绍F1大赛,我对当中一些人得的名字还是比较熟悉的。

新闻稿中说,自2001年雷诺车队成立以来,CD-Adapco就是雷诺车队的CFD软件提供商。RS26引擎R26底盘 的设计都和CD-Adapco的STAR-CD或STARE-CCM+软件有关。CD-Adapco是F1大赛中最成功的CFD软件提供商,因为,在2006赛季有得分的车队中,有70%的车队使用STAR-CD或者STAR-CCM+作气动或者引擎流动模拟。看来STAR-CD在汽车行业的优势是当之无愧的。关于F1大赛的更多信息可以参见搜狐体育的F1频道

星期日, 十一月 05, 2006

多面体网格的优点

多面体网格(polyhedral mesh)是最近CFD求解器开发的热点,它已经被使用在CD-Adapco的STAR-CCM+中,STAR-CD V4也支持多面体网格,Fluent V6.3也将支持多面体网格。其实在CFX中,由于采用基于节点的控制体积(Vertex-based control volume),用户提供给求解器的网站是正常的四面体或六面体网格,而在求解器内部会生成多面体网格,并且声称已经使用这个技术超过15年。下面的文章翻译自CD-Adapco上面的一个新闻发布(链接),作者是Milovan Peric和Stephen Ferguson。

自CD-Adapco推出STAR-CCM+求解器和自动多面体网格生成器以来, 我们经常被问到一个问题:为什么我们会采用多面体网格,而不是四面体网格? 在这篇文章里,我们会就这个问题提供答案和实例。

四面体单元是最简单的体积单元;他们的表面是平面碎片 (plane segment), 因此面和体的中心点很容易定义。四面体网格也是相对容易自动生成的。它们几乎是自动网格生成技术中的标准应用,主流CFD求解器都支持四面体网格。但是,在另一方面,四面体网格不能被拉伸太多,因此为了保证在边界层、长的管道、或小的缝隙的精确度, 需要使用比结构化网格(六面体)多的多的四面体网格。棱形网格层可以部分解决这个问题。

四面体控制体只有四个相邻网格,在使用标准近似(线性形函数)来计算单元中心的梯度时,会比较困难。一个困难是相邻网格节点的空间位置,因为它们可能会几乎处在一个平面上,使得几乎不可能计算垂直于该平面的梯度。当单元邻近边界时会有另外一个可能碰到的问题:即使只有一个面是边界面,其它三个面的分布可能会很不好;尤其是在邻近线或者角落时,一个单元可能只有两个或者一个相邻单元,这不光会降低精确度,而且会引起严重的数值问题。

为了在四面体网格上实现精确求解和好的收敛性,需要特别的离散技术和大量的单元。这些解决方案并不是最优化的;第一种方法使得程序更加复杂、更难扩展和维护,而第二个方法增加了内存和计算时间需求。多面体网格提供了和四面体网格一样的自动网格生成功能,并且能解决这些缺点。

多面体网格的主要优点是它有多个相邻单元(通常10数量级) ,和四面体网格相比,使用多面体网格能更好的近似梯度(使用线性形函数以及最近单元的信息)。即使靠近线和角,一个多面体单元也可能有几个相邻单元,从而可以合理的预测梯度和当地流动分布。更多的相邻单元意味着更多的内存使用量和更多的计算操作,这可以从更高的精度中得到补偿。下面将要讨论这个问题。

多面体网格也不如四面体网格多拉伸敏感。智能化的网格生成和优化技术也提供了无限的可能性:通过加入新的节点、单元或面,单元可以自动结合、分裂、或修改。实际上,网格质量在将来会有显著提高,并提高求解器效率和精确度。而且,求解器中在使用四面体网格时需要的特别处理将不再必要。在使用多面体网格时,要模拟当地网格细化、滑移网格界面、周期性边界条件等,只需要生成特别的多面体网格,对求解器来说它们都是一样的。

多面体网格在处理回流时更有优势。比如测试显示在模拟正方体空腔驱动流时,要得到指定的精度所需的多面体单元数量要比笛卡尔六面体(通常认为对长方体计算域,六面体是最理想的)少得多。这可以解释为:对六面体单元,有三个最佳流动方向可以达到最大精度(和三对平行面垂直的方向) ;对有十二个面的多面体单元,有六个最佳流动方向,并且有更多的相邻单元,因此可以用较少的单元数得到较高的精确度。对于不同网格类型特点和一个测试结果可以参考Peric (2004) 发表的文章。

很多实际应用的比较分析显示,要达到和四面体网格同样的精确度,使用多面体网格只需要四分之一的单元数、一半的内存、十分之一到五分之一的计算时间。另外,多面体网格有更好的收敛性,通常不需要调整默认的求解器参数,如图1至3所示。

图1. 模拟发动机水套内流动所用的多面体网格


图2. 最细的四面体网格(左图)和最细的多面体网格(右图)所预测的水套壁面压力分布

图3. 不同四面体网格和多面体网格所预测的水套入口出口压力差

这个测试模拟的是发动机水套;监测的是进口和出口的压力差。计算分别在6个多面体网格(单元数21872至593888)和6个四面体网格(单元数39587至2322106)上进行。在所有情形中近壁面处均使用了棱形网格层。图1是一个多面体网格。图2是在最细的多面体网格和最细的四面体网格上计算所得的压力分布。因为网格很小,这两个计算结果非常相似。为了分析结果对网格的依赖性,图3显示了在不同网格上计算得到的压差。该图显示这两种网格类型都可以收敛得到网格无关解。所有的情形使用了同样的离散和求解方法(对流项采用二阶迎风差分格式)。

重要的是,在多面体网格上得到的解总是比有相近网格数量的四面体网格要精确。实际上,从有65513单元的多面体网格上得到的解比从有393273个单元(大约6倍多)的四面体网格上得到的解还要稍微精确。多面体网格所需的计算时间只是得到同等精度的四面体网格的十分之一。到目前为止所进行的比较得到了相似的结果,这验证了CD-Adapco发展这一技术的决定是正确的。随着CD-Adapco的全自动多面体网格生成器的面世,多面体网格将会成为业界标准。

参考文献:

1 M. Peric: Flow simulation using control volumes of arbitrary polyhedral shape, ERCOFTAC Bulletin, No. 62, September 2004.

一周没有发言了

原因是参加了一个关于CAE的conference。

星期三, 十月 25, 2006

Abaqus 6.6学生版发布





达索系统(Dassault systems,全球网站中国网站)在本月10号发布了Abaqus
6.6的学生版(Abaqus全球网站中文网站新闻发布)。Abaqus学生版本是研究版本的教学版本的补充,对学习Abaqus进行有限元分析的学生和专业人士有很大的帮助。



在Abaqus网站上对产品的特色有详细介绍。主要的限制是用户不能求解和后处理超过1000个节点的模型。另外,用户不能编程扩展。Abaqus也原则上不提供技术支持。对软件的功能和模型没有特别的限制,也就是说用户可以接触练习所有Abaqus模型。价格是99美元一套。这样的价格应该是非常合理的。不知道ANSYS有没有提供这样的学生版本,如果有的话,我一定会买一套,但是在ANSYS的官方网站上并没有相关信息。看来ANSYS对教育市场还是比较保守的。

星期一, 十月 09, 2006

压力和应力的区别

今天有个朋友问我压力和应力的区别。压力(pressure)和应力(stress)有相同的单位帕斯卡(Pa, N/m^2),但是它们的区别是什么呢?


通常来说,固体内部只有应力,只有流体内部才既有应力又有压力。应力是物体对外力的方应,而压力和温度、压力一样是流体的一种状态。应力是矢量,有方向,而压力是标量没有方向。应力是由于材料的弹性或粘性而产生的,压力是由于流体内部分子的无规则热运动而产生的。


虽然在固体力学中,一些外力可以定义成压力。比如,当我们用手指摁桌面时,会给桌面施加一个压力。这里的压力是指和应力相对应物体内部的压力。


综合来说,应力是在外力作用下固体或流体由于材料的弹性或粘性对外力的反应。压力是由于流体内部的无规则分子热运动而产生的相互作用力,和材料的粘性无关。


压力的定义是相对简单的(wikipediaNASA Glenn Research Center )。因为压力是由于无规则分子热运动产生的,因此各向同性,是个标量。而应力的定义就要复杂得多。


1. 应力是表面力


当一个物体受到外力作用时,物体内部会相应产生内力来平衡外力。而内力的分布和该物体的强度以及断裂等直接相关。内应力(stress)就是是描述物体内部力的分布的矢量。值得注意的是,应力是表面力(surface force,单位面积上的力),而不是体积力(volume force,单位体积上的力,比如重力)。如果我们把一个材料划分成无限小的连续微粒,微粒之间的相互作用都是表面力(因为微粒在面上接触)。对于单个微粒而言,在静力学中,由于受力平衡,体积力为零;在动力学中,体积力正比于加速度。显然,为了表示材料的强度,应该使用表面力。而体积力可以应用在动力学上,比如在流体动力学中,动量方程的左边是体积力。


2. 应力是对外力的反应



如上图所示,一个外力(traction force)作用在任意物体表面上。要完整表示这个外力,必须既要描述这个力(大小和方向)和该力所作用的表面(面积和方向)。比如,只给定外力的大小和方向,而没有给定作用表面的面积和方向,该力的作用效果是不确定的(或者说定义不完整)。


作用在一个无限小表面上的外力,T,可以表示为:



外力T是一个二维张量,表示无限小表面上表面力的大小,单位是帕斯卡。



应力是对外力的反应。我们可以用同样的方法来定义内应力。假设如上图所示,一个物体被任意曲面剖开。在该剖面上任意无限小表面上的应力可以表示为



因此,应力可以理解为在物体内部平面上的内力的分布。要描述应力,必须要先定义一个平面。


3. 物体内一点的应力




应力是物体内部平面上的内力分布,通常用二阶张量来表示物体内一点(无限小的连续微粒)的应力。


在笛卡尔坐标系中,物体内一个无限小微粒的应力分布可以用如上图所示的无限小正方体来表示(在笛卡尔坐标系中,用正方体可以得到最简洁的形式)。在正方体的六个表面中,每个表面有三个相互正交的应力分量(一个正应力分量,,两个切应力分量,)。显然,和其他矢量一样,各分量的大小是和坐标系的选取有关的。




在该正方体表面上,实际有6×3=18个分量,因此这种用正方体表示应力的方法容易引起混淆。当这个无限小的正方体收缩成一个点时,该点的应力可以用下图来表示。通过该点有三个正交平面,每个平面上有三个分量,共有3×3=9个分量。这样,我们就得到了小正方体图中的矩阵:


在该矩阵中,代表正应力(垂直于作用平面),代表切应力(平行于作用平面),第一个下标代表应力作用平面的法线方向,第二个下标代表应力分量的方向。比如代表作用在x平面,y方向的切应力。



值得注意的是,上面的矩阵总是相对于对角线对称的。也就是说,。这是因为微粒表面上的力矩是平衡的,否则该微粒将会有无限大的旋转加速度。如上图所示的二维情况,如果,该微粒表面上的力相对于A点的力矩为(假设正方形边长为1,力矩的定义为)。要满足力矩平衡,必须有


因此,应力可以由六个独立的分量表示:



至此,我们已经得到了最常见的应力的定义。